R语言中复杂目标函数定义与大规模列数下MIP模型不可行问题咨询

Hey there, let's work through your questions step by step:

1. Can we define complex objective functions in R optimizers? #

Absolutely! For mixed-integer programming (MIP) tasks like yours, R has great tools—think ompr, ROI, or lpSolve—that let you translate logical objectives (like counting how many columns hit a threshold) into math the solver can handle. The trick is turning your "countif" logic into linear constraints using binary variables, which is exactly what you're trying to do.

2. R tools for maximizing the count of columns where sum exceeds target #

Yep, the tools you're already using (ompr + ROI) are perfect for this. The standard playbook here is to create a binary variable A[j] for each column j: set it to 1 if the sum of your selected rows for that column meets or beats the target, 0 otherwise. Then your objective is just maximizing the sum of all A[j]—that's exactly your "count the number of columns that meet the target" goal. Even with 10,000 columns, this approach works (you'll just need to pay attention to solver speed and model efficiency as you scale).

3. Troubleshooting your infeasible model at scale #

Now let's dig into why your model works for small column counts but goes infeasible at 150 columns, and fix it up:

Why your current model fails #

Looking at your constraint logic, there's a critical mix-up in what A[j] is supposed to represent. Let's break down what happens when A[j] = 0 (meaning you don't want that column to meet the target):

• The constraints s[j] <= 9999999*A[j] and s[j] >= A[j] force s[j] = 0.
• That turns your main column constraint into sum(vals[i,j]*x[i]) >= targets[j].

Wait a second—that means every column you mark as not meeting the target is actually being forced to meet the target! That's a huge contradiction. When you scale up to 150 columns, it's impossible for most columns to hit their target with only 5 selected rows, so the solver can't find a solution that satisfies all these conflicting rules—hence the "infeasible" result.

A corrected, more efficient model #

Let's rewrite the model to correctly capture your goal: A[j] = 1 if the column sum meets the target, 0 otherwise. Here's how to do it right:

1. Keep your core variables: x[i] (binary, 1 if row i is selected) and A[j] (binary, 1 if column j hits the target).
2. Keep the row selection constraint: sum(x[i]) = 5 (or 4 in your small example).
3. Link A[j] to the column sum with Big M:
   • First, calculate a reasonable Big M value—it should be larger than the maximum possible sum any column could have (sum of all rows in that column, since you're selecting up to 5 rows). This ensures the constraint doesn't artificially restrict valid solutions.
   • Add two linking constraints per column:
     • One that enforces: if A[j] = 1, the column sum must be >= target.
     • Another that enforces: if A[j] = 0, the column sum must be < target.
4. Set the right objective: Maximize the sum of A[j]—this directly counts how many columns meet your target.

Here's the revised code for your small example:

nr <- 10 # number of rows
nt <- 15 # number of target columns
vals <- matrix(sample.int(nr*nt, nr*nt), nrow=nr, ncol=nt)
targets <- vector(length=nt)
targets[1:nt] <- 4*mean(vals)

# Calculate a sensible Big M: max possible sum for any column (sum of all rows in the column)
M <- max(colSums(vals))

model <- MIPModel() %>% 
  add_variable(x[i], i = 1:nr, type = "binary") %>% 
  add_constraint(sum_expr(x[i], i = 1:nr) == 4) %>%  # Select 4 rows
  add_variable(A[j], j = 1:nt, type = "binary") %>% 
  # If A[j] = 1, column sum >= target; if A[j] = 0, no restriction here
  add_constraint(sum_expr(vals[i,j]*x[i], i = 1:nr) >= targets[j] - M*(1 - A[j]), j = 1:nt) %>% 
  # If A[j] = 0, column sum < target; if A[j] = 1, no restriction here
  add_constraint(sum_expr(vals[i,j]*x[i], i = 1:nr) <= targets[j] - 1e-6 + M*A[j], j = 1:nt) %>% 
  set_objective(sum_expr(A[j], j = 1:nt), "max")  # Maximize count of columns meeting target

# Solve the model
model <- solve_model(model, with_ROI(solver = "glpk"))

Tips for scaling to 10,000 columns #

• Upgrade your solver: GLPK is great for small tests, but for 10k columns, you'll want a faster solver. Commercial options like Gurobi or CPLEX (both have R interfaces) are top-tier, but open-source alternatives like CBC (via ROI.plugin.cbc) are way better than GLPK for large models.
• Tune solver settings: Enable presolve, adjust cut generation, or set a reasonable time limit to get feasible solutions faster—large MIPs can take time, so these tweaks help.
• Cut redundant variables: Your original model had s[j] which isn't needed in the revised approach—removing extra variables shrinks the model size and speeds up solving.

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