Validate solver solution count before processing

[josalhor]

Add validation for the number of solutions returned by the solver.

Should fix #876

This PR has not been tested (and I'm actually a little bit unsure why there is the try catch in the first place).
But the solution probably goes in this direction.

[CLAassistant]

All committers have signed the CLA.

[pchtsp]

thanks, can you provide a unit test with a toy example that currently fails and should pass?

[josalhor]

I think this should work:

import pytest
import pulp as lp
from pulp import LpProblem, LpMinimize, LpVariable, lpSum, LpBinary, LpInteger, LpStatus


TL = 4.0
TOL = 1e-6
BIG = 99999


N = ['n0','n1','n2','n3','n4','n5','n6']      
S = ['n4','n5','n6']                           
H = 10                                         
M = ['m0','m1']                                
K = ['k0','k1']                                
Q = ['q0','q1','q2']                           
W = {'n2':{'q0': 1.1,'q1':1.3,'q2':0.9}, 'n3':{'q0':1.2,'q1':1.2,'q2':0.9}}
R = {'m0':'n0','m1':'n1'}
C = {'m0':{'t0':100,'t1':120,'t2':260}, 'm1':{'t0':100,'t1':120,'t2':260}}
CAP = {'m0':100,'m1':100}
D = {'n4':{'q0':145,'q1':260,'q2':260}, 'n5':{'q0':0,'q1':370,'q2':75}, 'n6':{'q0':0,'q1':295,'q2':150}}
AL = 0.05

E = {
 ('n0','n2'):5, ('n2','n4'):6, ('n4','n5'):4, ('n0','n3'):6,
 ('n5','n6'):5, ('n6','n2'):12, ('n4','n6'):8, ('n3','n4'):6,
 ('n0','n4'):10, ('n0','n5'):15, ('n0','n6'):20, ('n3','n5'):4, ('n3','n6'):4,
 ('n2','n5'):9, ('n2','n6'):12, ('n3','n2'):6,
 ('n1','n2'):7, ('n1','n3'):5, ('n1','n4'):10, ('n1','n5'):8, ('n1','n6'):14, ('n1','n0'):10,
}
T = E.copy()
for (i,j),v in list(E.items()):
    E[j,i]=v; T[j,i]=v
for i in N:
    E[i,i]=999; T[i,i]=0
for _,r in R.items():
    E[r,r]=0

def build():
    P = LpProblem("p", LpMinimize)
    x = LpVariable.dicts("x", ((m,i,j,t) for m in M for i in N for j in N for t in range(H-1)), 0, 1, LpBinary)
    y = LpVariable.dicts("y", ((m,i,t) for m in M for i in N for t in range(H)), 0, 1, LpBinary)
    z = LpVariable.dicts("z", ((m,t,q,c) for m in M for t in range(H) for q in Q for c in C[m]), 0, None, LpInteger)
    u = LpVariable.dicts("u", ((m,t,q,c,n) for m in M for t in range(H) for q in Q for c in C[m] for n in ['n2','n3']), 0, None, LpInteger)
    d = LpVariable.dicts("d", ((m,s,q,c,t) for m in M for s in S for q in Q for c in C[m] for t in range(H)), 0, None, LpInteger)

    P += lpSum(E[i,j]*x[m,i,j,t] for m,i,j,t in x) + lpSum(u[m,t,q,c,n]*W[n][q] for m in M for t in range(H-1) for q in Q for c in C[m] for n in ['n2','n3'])

    for m,r in R.items():
        P += y[m,r,0] == 1
    for m in M:
        for q in Q:
            for c in C[m]:
                P += z[m,0,q,c] == 0
                for n in ['n2','n3']:
                    P += u[m,0,q,c,n] == 0

    for m in M:
        P += lpSum(T[i,j]*x[m,i,j,t] for i in N for j in N for t in range(H-1)) + \
             lpSum(u[m,t,q,c,n]*AL for t in range(H-1) for q in Q for c in C[m] for n in ['n2','n3']) + \
             lpSum(d[m,s,q,c,t]*AL for s in S for q in Q for c in C[m] for t in range(H-1)) <= CAP[m]

    for m,r in R.items():
        P += lpSum(x[m,i,r,t] for t in range(H-1) for i in N if i!=r) == 1
        P += lpSum(x[m,r,i,t] for t in range(H-1) for i in N if i!=r) == 1
        P += y[m,r,H-1] == 1

    for m in M:
        for t in range(H-1):
            for j in N:
                P += y[m,j,t+1] == lpSum(x[m,i,j,t] for i in N)
            for i in N:
                P += lpSum(x[m,i,j,t] for j in N) == y[m,i,t]
        for t in range(H):
            P += lpSum(y[m,i,t] for i in N) == 1

    for m in M:
        for t in range(1,H):
            for q in Q:
                for c in C[m]:
                    P += z[m,t,q,c] == z[m,t-1,q,c] - lpSum(d[m,s,q,c,t] for s in S) + lpSum(u[m,t,q,c,n] for n in ['n2','n3'])
                    P += z[m,t,q,c] <= C[m][c]

    for m in M:
        for t in range(1,H):
            for q in Q:
                for c in C[m]:
                    for n in ['n2','n3']:
                        P += u[m,t,q,c,n] <= y[m,n,t]*C[m][c]

    for m,r in R.items():
        for t in range(1,H):
            P += BIG*(1 - y[m,r,t]) >= lpSum(z[m,t-1,q,c] for q in Q for c in C[m])
            for n in ['n2','n3']:
                P += BIG*(1 - y[m,n,t]) >= lpSum(z[m,t-1,q,c] for q in Q for c in C[m])

    for m in M:
        for s in S:
            for t in range(H):
                for q in Q:
                    P += lpSum(d[m,s,q,c,t] for c in C[m]) <= y[m,s,t]*D[s][q]

    for s in S:
        for q in Q:
            P += lpSum(d[m,s,q,c,t] for m in M for c in C[m] for t in range(H)) == D[s][q]

    for m in M:
        for t in range(1,H):
            for q in Q:
                for c in C[m]:
                    P += lpSum(d[m,s,q,c,t] for s in S) <= z[m,t-1,q,c]

    return P

def _claimed(prob):
    s = LpStatus.get(prob.status, str(prob.status))
    return s in {"Optimal"} or "Feasible" in s

def _check_integral_and_constraints(prob):
    
    for v in prob.variables():
        val = v.varValue
        if val is None:  
            continue
        if v.cat == "Binary":
            if abs(val - round(val)) > TOL or round(val) not in (0,1):
                return False
        elif v.cat == "Integer":
            if abs(val - round(val)) > TOL:
                return False
    
    for c in prob.constraints.values():
        lhs = c.value()
        rhs = c.constant
        if lhs is None:
            return False
        if c.sense == lp.LpConstraintEQ and abs(lhs - rhs) > TOL: return False
        if c.sense == lp.LpConstraintLE and lhs - rhs > TOL: return False
        if c.sense == lp.LpConstraintGE and rhs - lhs > TOL: return False
    return True

def _solve_scip(prob):
    prob.solve(lp.SCIP_PY(msg=True, options=["limits/time", TL]))

@pytest.mark.parametrize("tlim", [TL])
def test_scip_timelimit_behavior(tlim):
    p = build()
    _solve_scip(p)
    s = LpStatus.get(p.status, str(p.status))

    if _claimed(p):
        ok = _check_integral_and_constraints(p)
        assert ok, "solver claimed a model but it does not satisfy the model"
    else:
        assert any(x in s for x in ["Infeasible", "Not Solved", "Undefined", "No Solution", "Time"]), \
               f"unexpected status without model: {s}"

if __name__ == "__main__":
    p = build()
    _solve_scip(p)
    
    raise AssertionError("run with pytest; this script is designed to fail under __main__")

The problem really does not matter, what matters is that the solver stops in the time limit before it has found a solution.