Thursday, November 18, 2021

Solve status: not as easy as it looks

In [1] a question was asked about a solve status in CVXR when we hit a time limit. Hopefully, we obtain the best solution found so far. But do we also get some indication about the status of this solution? CVXR only has the following status codes [2]:

status The status of the solution. Can be "optimal", "optimal_inaccurate", "infeasible", "infeasible_inaccurate", "unbounded", "unbounded_inaccurate", or "solver_error".

That looks wholly inadequate to properly describe this situation.

I tried a few scenarios for a MIP model.

  1. Optimal. There is no problem with this.  I get:

    > result$status
    [1]   "optimal"


  2. Time limit, no integer solution found.

    > result$status
    [1]   "unbounded"


    This is somewhat of a surprise. Anyway, it is wrong. However, there is no good status that would describe this situation.

  3. Time limit, integer solution found.

    > result$status
    [1] "optimal"            "optimal_inaccurate"
    Warning messages:
    1: In if (solution[[STATUS]] %in% SOLUTION_PRESENT) { :
      the condition has length > 1 and only the first element will be used
    2: In if (status %in% SOLUTION_PRESENT) { :
      the condition has length > 1 and only the first element will be used
    3: In if (!(solution@status %in% SOLUTION_PRESENT)) return(Solution(solution@status,  :
      the condition has length > 1 and only the first element will be used
    4: In if (solution@status %in% SOLUTION_PRESENT) { :
      the condition has length > 1 and only the first element will be used


    This is problematic. I see two status codes and none of them really covers this case. There are also a number of warning messages. This just looks like a bug.

Obviously, we see some problems here. Some of them may be just bugs, but there is also a conceptual problem: CVXR status codes have a way too simplistic view of how a solver may terminate. GAMS almost takes the opposite view. It has two return codes:


Model StatusSolver Status
1OPTIMAL1NORMAL COMPLETION
2LOCALLY OPTIMAL2ITERATION INTERRUPT
3UNBOUNDED3RESOURCE INTERRUPT
4INFEASIBLE4TERMINATED BY SOLVER
5LOCALLY INFEASIBLE5EVALUATION INTERRUPT
6INTERMEDIATE INFEASIBLE6CAPABILITY PROBLEMS
7FEASIBLE SOLUTION7LICENSING PROBLEMS
8INTEGER SOLUTION8USER INTERRUPT
9INTERMEDIATE NON-INTEGER9ERROR SETUP FAILURE
10INTEGER INFEASIBLE10ERROR SOLVER FAILURE
12LIC PROBLEM - NO SOLUTION11ERROR INTERNAL SOLVER FAILURE
13ERROR NO SOLUTION12SOLVE PROCESSING SKIPPED
14NO SOLUTION RETURNED13ERROR SYSTEM FAILURE
15SOLVED UNIQUE
16SOLVED
17SOLVED SINGULAR
18UNBOUNDED - NO SOLUTION
19INFEASIBLE - NO SOLUTION

The total number of combinations is quite large. Of course, this still does not cover all possible return codes. A famously confusing return code is "Model is infeasible or unbounded" (Gurobi, without a solution). This cannot be mapped perfectly to any of these codes.

One of the important things to learn from these codes is that there are different forms of "infeasible". The solver may not have found a feasible solution yet (this is called "intermediate infeasible" in the above table). Or it may be proven to be infeasible. Or a local nonlinear solver may not be able to find a feasible solution, but it can not prove there are none ("locally infeasible").

What can we expect in a MIP solver?


When solving a MIP, the model goes through different phases. Each can fail or be interrupted. In general, it can be assumed the best available solution so far will be returned. Solvers don't follow this rule as religious as I would like. E.g. when there is no integer solution yet, a solver could return the relaxed solution. Solvers often do not do this. 

References


  1. CVXR in R: What happens when no solution is found within the time limit?, https://stackoverflow.com/questions/69982620/cvxr-in-r-what-happens-when-no-solution-is-found-within-the-time-limit
  2. Package CVXR, https://cran.r-project.org/web/packages/CVXR/CVXR.pdf


Appendix: model code


This is the code I used to test this:

library(CVXR)

installed_solvers()

set.seed(1234)

m <- 100
n <- 50

A <- matrix(runif(m*n,min=0,max=10),nrow=m,ncol=n)
b <- runif(m,min=10,max=20)


x <- Variable(n)
r <- Variable(m)
d <- Variable(n,boolean=T)
k <- Parameter()
obj <- Minimize(sum(abs(r)))
constr <-list(A %*% x - b == r, 
              x >= -d*1000, 
              x <= d*1000, 
              sum(d) == k)
prob <- Problem(obj,constr)

# find optimal solution
value(k) <- 2 # easy problem
result <- solve(prob, solver="GLPK", verbose=T)
result$status

# time limit before we have an integer solution
value(k) <- 10 # difficult problem
result <- solve(prob, solver="GLPK", verbose=T, tm_limit = 1)
result$status

# find integer solution (non-optimal)
result <- solve(prob, solver="GLPK", verbose=T, tm_limit = 1000)
result$status

1 comment:

  1. This is unfortunately a known issue in the CVX family (had the same problem with CVXPY); not sure why that is. Maybe a pull request is due.

    ReplyDelete