Termination condition is optimal but number of solutions is 0

[tristan-bifrost]

I am using pyomo (Pyomo 5.7.3 (CPython 3.9.6 on Windows 10)) with cbc (CBC 2.10.5) to solve a MILP. Using the following command

results = opt.solve(self.model, tee=False, logfile="model.log",options_string=solver_parameters)
print(results)

I get the following

Problem:
- Lower bound: -inf
  Upper bound: inf
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: unknown
Solver:
- Status: ok
  Message: CBC 2.10.5
  Termination condition: optimal
  Id: 0
  Error rc: 0
  Time: 0.012964248657226562
Solution:
- number of solutions: 0
  number of solutions displayed: 0

Status if ok and termination is optimal, however the number of solution if 0, which is confirmed when trying to access the objective value. I've used different solvers with pyomo, as well as stand alone cbc, and, and in any case the solver returned a status of infeasible.

I'm also surprised that number of constraints and number of variables are 0.

[jsiirola]

Please include a minimal example that reproduces the problem. It looks like the SOL file parser is not picking up the "infeasible" message from that version of CBC.

Also note that Pyomo 5.7.3 is more than a year old. Please consider upgrading and testing against the current release.

[tristan-bifrost]

Updating to pyomo 6.2 doesn't solve the issue. I've tried using several versions of cbc, the problem remains.

I initially posted the question on OR stackexchange : https://or.stackexchange.com/questions/7810/problem-is-infeasible-with-gurobi-feasible-with-cbc-but-cant-access-objective The question contains lp and mps files. It is difficult to provide the complete code to reproduce de problem because it is part of a large project.

[blnicho]

The title of this issue is a duplicate of #556.

By default, Pyomo loads the solution returned by the solver directly back into the model (it is more efficient than generating the representation used by the results object), so the results object (correctly) reports that it is holding no solutions. You can see the model status with:

model.display()

Originally posted by @jsiirola in https://github.com/Pyomo/pyomo/issues/556#issuecomment-395742933

Your comment about the number of constraints and variables being 0 is also duplicated in #1033

As for issues with different versions of CBC, there are several known issues and feature requests summarized in #1926.

I'm going to close this as I think all of your bug reports are captured in the issues linked above.

[tristan-bifrost]

I don't see how #556 relates to my problem, which is not about printing information, but about pyomo returning the wrong status and termination condition.

I don't see how any request of #1926 relates to the fact that pyomo is returning the wrong status and termination condition with cbc.

I don't understand how you answered the issue and why you closed it.

Again, the issue is pyomo returning ok and optimal when it should return infeasible.

[blnicho]

If you solve the model with the tee=True option, what is the termination status shown in the solver log? If there is an inconsistency in the termination status reported in the solver log and the one reported by Pyomo we certainly want to know about it but you haven't included a solver log in your report and the title of this issue is not reporting that as the problem.

If you read through issues #556 and #1033 they describe why printing the results object doesn't really give you the information you might expect it to and is not a recommended way for querying the result. Using model.display or check_optimal_termination(results) are both better ways to check the solver results loaded into Pyomo.

If I have misinterpreted your issue report I apologize, feel free to reopen it. Including a snippet of the solver log showing the termination status would help us pinpoint the issue.

[michaelbynum]

I put together a small infeasible problem (see below), and Pyomo correctly reports the problem as infeasible. However, this could certainly depend on the versions of Pyomo and CBC. I was using Pyomo 6.2 and CBC 2.10.7. I agree with @blnicho - more information is probably needed to debug this. It sounds like an edge case that may be difficult to replicate.

In [1]: import pyomo.environ as pe

In [2]: m = pe.ConcreteModel()

In [3]: m.x = pe.Var(bounds=(3, 5))

In [4]: m.y = pe.Var()

In [5]: m.obj = pe.Objective(expr=m.y)

In [6]: m.c1 = pe.Constraint(expr=m.y >= m.x)

In [7]: m.c2 = pe.Constraint(expr=m.y <= -m.x)

In [8]: opt = pe.SolverFactory("cbc")

In [9]: res = opt.solve(m, tee=True)
Welcome to the CBC MILP Solver
Version: 2.10.7
Build Date: Jan 26 2022

Option for printingOptions changed from normal to all
Presolve determined that the problem was infeasible with tolerance of 1e-08
Presolved model looks infeasible - will use unpresolved

Problem has 3 rows, 3 columns (1 with objective) and 5 elements
There are 1 singletons with no objective
Column breakdown:
1 of type 0.0->inf, 0 of type 0.0->up, 0 of type lo->inf,
1 of type lo->up, 1 of type free, 0 of type fixed,
0 of type -inf->0.0, 0 of type -inf->up, 0 of type 0.0->1.0
Row breakdown:
0 of type E 0.0, 1 of type E 1.0, 0 of type E -1.0,
0 of type E other, 0 of type G 0.0, 0 of type G 1.0,
0 of type G other, 2 of type L 0.0, 0 of type L 1.0,
0 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other,
0 of type Free
Presolve determined that the problem was infeasible with tolerance of 1e-08
Analysis indicates model infeasible or unbounded
0  Obj 0 Primal inf 6.9999997 (3) Dual inf 0.0099999 (1) w.o. free dual inf (0)
3  Obj 0 Primal inf 2.9999999 (1) Dual inf 0.4999999 (1)
4  Obj 3 Primal inf 5.9999999 (1)
Primal infeasible - objective value 3
PrimalInfeasible objective 3 - 4 iterations time 0.002

Result - Linear relaxation infeasible

Enumerated nodes:           0
Total iterations:           0
Time (CPU seconds):         0.00
Time (Wallclock Seconds):   0.00

Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

WARNING: Loading a SolverResults object with a warning status into
    model.name="unknown";
      - termination condition: infeasible
      - message from solver: <undefined>

In [10]: res.solver.termination_condition
Out[10]: <TerminationCondition.infeasible: 'infeasible'>

[jsiirola]

I started diving into the CBC source. There appears to be a set of infeasible / not optimal messages that can be emitted that I do not think are caught by our parser (emitted here) [and I think these are what are being emitted by this example]. That said, it is really difficult to see from the code how those messages get triggered, so I am not sure how to generate a Pyomo test case that exercises this particular branch of the CBC code.

[michaelbynum]

Interesting...

[jsiirola]

I converted the LP file from the SO ticket into an equivalent (flat) Pyomo model (issue-2290.py.txt). When sent to CBC, the interface correctly identifies the model as Infeasible (in the presolve).

At this point, I am not sure there is much more we can do without a Pyomo model that reproduces the problem.

[jessecurry]

I seem to be running into the issue described here, though I'm probably doing something stupid as this is my first Pyomo project.

Here is a model that reproduces the problem.

import os
import numpy as np
import pandas as pd
import pyomo.environ as pe

#####
# Constants
days = ["monday", "tuesday", "wednesday", "thursday", "friday", "saturday", "sunday"]
activities = ["hard", "moderate", "active-recovery", "low-impact"]
activityCost = np.array([0.8, 0.85, 0.9, 0.75])
activityCapCost = np.array([0.5, 0.25, 0.1, 0.75])  # Will be added to 1 and applied

# Cost matrices
activityCosts = np.matrix(np.tile(activityCost.transpose(), (len(days), 1))).transpose()
activityCapCosts = np.matrix(
    np.tile(activityCapCost.transpose(), (len(days), 1))
).transpose()
activityAffinities = np.matrix(np.ones((len(activities), len(days))))
# activityAffinities[1, 5] = 0.015
# activityAffinities[1, 6] = 0.015

#####
# Inputs
previousDay = [0, 0, 0, 0]
activityReq = [2, 2, 0, 0]
remainingDays = 7

#####
# Model
model = pe.ConcreteModel("activity Scheduler")
model.activities = pe.Var(
    ((day, activity) for day in days for activity in activities),
    within=pe.Binary,
    initialize=0,
    bounds=(0, 1),
)

#####
# Objective
def objective_rule(model, verbose=False):
    # Build activity matrix and capacity matrix
    activityMtx = np.matrix(np.zeros((len(activities), len(days))))
    capMtx = np.matrix(np.zeros((len(activities), len(days))))
    capMtx[:, 0] = np.reshape(previousDay, (len(activities), 1))
    for i, day in enumerate(days):
        for j, activity in enumerate(activities):
            activityMtx[j, i] = pe.value(model.activities[day, activity])
            if i < 6:
                capMtx[j, i + 1] = pe.value(model.activities[day, activity])

    # Calculate basic activity costs
    cost_matrix = np.multiply(activityMtx, activityCosts)

    # Apply affinity credit
    cost_matrix = np.multiply(cost_matrix, activityAffinities)

    # Add penalty costs
    penalty_matrix = np.multiply(capMtx, activityCapCosts)
    penalty_vector = penalty_matrix.sum(axis=0)
    penalty_vector = np.add(penalty_vector, 1.0)
    cost_matrix = np.multiply(cost_matrix, penalty_vector)

    return cost_matrix.sum()

model.obj = pe.Objective(
    rule=objective_rule, sense=pe.minimize, name="activity Schedule"
)

#####
# Constraints
model.constraints = pe.ConstraintList()

# No more than one activity per day
for day in days:
    model.constraints.add(
        sum(model.activities[day, activity] for activity in activities) <= 1
    )

# Do not exceed the activity requirements
for i, activity in enumerate(activities):
    model.constraints.add(
        sum(model.activities[day, activity] for day in days) <= activityReq[i]
    )

# Assign the maximum allowable activities
model.constraints.add(
    sum(model.activities[day, activity] for day in days for activity in activities)
    == min(remainingDays, sum(activityReq))
)

#####
# Solution
solver_manager = pe.SolverManagerFactory("serial")
solver = pe.SolverFactory("cbc")
# solver.options = {"maxIt": 268435456}

results = solver_manager.solve(model, opt=solver, tee=True)

# model.display()
print(f"Is optimal: {pe.check_optimal_termination(results)}")
# print(results)

activities_table = pd.DataFrame(
    np.zeros((len(activities), len(days))), index=activities, columns=days
)

for i, day in enumerate(days):
    for j, activity in enumerate(activities):
        activities_table.iloc[j, i] = pe.value(model.activities[day, activity])

def map_activity_type(day):
    if day[0] == 1:
        return "hard"
    elif day[1] == 1:
        return "moderate"
    elif day[2] == 1:
        return "active-recovery"
    elif day[3] == 1:
        return "low-impact"
    else:
        return None

print(activities_table.apply(map_activity_type, axis=0))

The basic idea here is that we're trying to optimally schedule a number of activities. Doing an activity two days in a row incurs a penalty. We want to bias towards scheduling harder activities, but have a higher penalty for activities after a hard day.

[jessecurry]

This is the result when I run this locally:

Welcome to the CBC MILP Solver 
Version: 2.10.7 
Build Date: Jul  6 2023 

command line - /opt/homebrew/opt/cbc/bin/cbc -printingOptions all -import /var/folders/sk/p45_mh7n1d7czj88jxqm_s3c0000gn/T/tmpuyev5vuh.pyomo.lp -stat=1 -solve -solu /var/folders/sk/p45_mh7n1d7czj88jxqm_s3c0000gn/T/tmpuyev5vuh.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 10 (-2) rows, 14 (-15) columns and 42 (-42) elements
Statistics for presolved model
Original problem has 28 integers (28 of which binary)
Presolved problem has 14 integers (14 of which binary)
==== 14 zero objective 1 different
14 variables have objective of 0
==== absolute objective values 1 different
14 variables have objective of 0
==== for integers 14 zero objective 1 different
14 variables have objective of 0
==== for integers absolute objective values 1 different
14 variables have objective of 0
===== end objective counts

Problem has 10 rows, 14 columns (0 with objective) and 42 elements
Column breakdown:
0 of type 0.0->inf, 0 of type 0.0->up, 0 of type lo->inf, 
0 of type lo->up, 0 of type free, 0 of type fixed, 
0 of type -inf->0.0, 0 of type -inf->up, 14 of type 0.0->1.0 
Row breakdown:
0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, 
1 of type E other, 0 of type G 0.0, 0 of type G 1.0, 
0 of type G other, 0 of type L 0.0, 7 of type L 1.0, 
2 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, 
0 of type Free 
Continuous objective value is 0 - 0.00 seconds
Cgl0002I 14 variables fixed
Cgl0004I processed model has 10 rows, 14 columns (14 integer (14 of which binary)) and 42 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 0
Cbc0038I Before mini branch and bound, 14 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 0 - took 0.00 seconds
Cbc0012I Integer solution of 0 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 0, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 0 to 0
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                0.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.01

Total time (CPU seconds):       0.00   (Wallclock seconds):       0.01

Is optimal: True
monday           None
tuesday      moderate
wednesday    moderate
thursday         hard
friday           None
saturday         hard
sunday           None

[blnicho]

@jessecurry if you would like to report a bug please open a new issue following our issue template rather than commenting on an old closed issue. If you think your bug report is related to a closed issue then you can easily link to it by including the issue number in your bug report.

[jessecurry]

@blnicho will do, wasn't sure how you all treat it and was responding to @jsiirola's comment:

At this point, I am not sure there is much more we can do without a Pyomo model that reproduces the problem.