在R中求解含全局函数约束的二进制变量约束优化问题

Yep, you hit a key limitation of MILP frameworks like ompr here—they only work with linear (or linearizable) constraints, and a global function f(x) can't be plugged directly into an add_constraint() call. Let's break down your solutions based on the nature of your f(x):

If your f(x) can be rewritten as a set of linear or mixed-integer linear constraints, this is the best path—you’ll get an exact solution using mature MILP solvers like GLPK or Gurobi.

For example, suppose f(x) is the sum of squared weights of selected elements: f(x) = sum((w[i] * x[i])^2). This looks non-linear, but since x[i] is binary (0 or 1), we can linearize it with auxiliary variables:

• Create a continuous variable z[i] where z[i] = w[i]^2 * x[i] (when x[i]=1, z[i] equals the squared weight; when x[i]=0, z[i]=0).
• Then your constraint becomes sum(z[i]) <= k, which is fully linear.

Here’s the adjusted ompr code:

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

n <- 10
w <- 1:n
k <- 60

model <- MILPModel() %>%
  add_variable(x[i], i = 1:n, type = "binary") %>%
  add_variable(z[i], i = 1:n, type = "continuous", lb = 0) %>%
  # Link z[i] to x[i] linearly
  add_constraint(z[i] == w[i]^2 * x[i], i = 1:n) %>%
  set_objective(sum_expr(x[i], i = 1:n), sense = "min") %>%
  add_constraint(sum_expr(z[i], i = 1:n) <= k)

result <- solve_model(model, with_ROI(solver = "glpk"))
get_solution(result, x[i])

Even for logical constraints (like "at least 3 consecutive 1s in x"), you can use standard MILP tricks to translate them into linear inequalities.

If f(x) is a black-box function, non-linear, or too complex to linearize, heuristic methods like genetic algorithms or simulated annealing can find good approximate solutions without needing linear constraints.

Let’s use the GA package for a genetic algorithm example:

library(GA)

n <- 10
k <- 60

# Define your global f(x) (example: product of selected indices)
f_x <- function(x) {
  selected <- which(x == 1)
  if (length(selected) == 0) return(0)
  prod(selected)
}

# Fitness function: maximize -sum(x) (since GA maximizes fitness; we want minimal sum(x))
fitness <- function(x) {
  if (f_x(x) <= k) {
    return(-sum(x))  # Valid solutions get higher fitness if sum(x) is smaller
  } else {
    return(-Inf)  # Invalidate solutions that violate the constraint
  }
}

# Run the genetic algorithm
ga_result <- ga(type = "binary", fitness = fitness, nBits = n, 
                popSize = 50, maxiter = 100, run = 10)

# Extract the best solution
best_x <- ga_result@solution[1,]
cat("Best x vector:", best_x, "\n")
cat("sum(x) =", sum(best_x), "\n")
cat("f(x) =", f_x(best_x), "\n")

This approach works for any computable f(x), but note it doesn’t guarantee a global optimal solution—just a strong approximation.

If n is small (say, n ≤ 20), you can brute-force all possible binary combinations, filter those that satisfy f(x) ≤ k, and pick the one with the smallest sum(x).

Here’s how to do this with purrr and expand.grid:

library(purrr)

n <- 10
k <- 60

# Generate every possible binary vector of length n
all_x <- expand.grid(rep(list(0:1), n))

# Define your f(x) (same product example)
f_x <- function(x) {
  selected <- which(x == 1)
  if (length(selected) == 0) return(0)
  prod(selected)
}

# Filter valid solutions and find the one with minimal sum(x)
valid_solutions <- all_x %>%
  mutate(
    f_val = pmap_dbl(., ~f_x(c(...))),
    sum_x = rowSums(.)
  ) %>%
  filter(f_val <= k)

best_solution <- valid_solutions %>%
  arrange(sum_x) %>%
  slice(1)

cat("Optimal sum(x):", best_solution$sum_x, "\n")
cat("Corresponding x vector:", unlist(best_solution[,1:n]), "\n")
cat("f(x) value:", best_solution$f_val, "\n")

This gives you a guaranteed global optimal solution, but becomes infeasible once n gets large (2^20 is over a million combinations, 2^30 is 1 billion—way too slow).

内容的提问来源于stack exchange，提问作者Jackk