Pyomo + GLPK does not return a certificate of optimality for simple LP?

[makansij]

Summary

I’m using Pyomo + GLPK to solve an LP involving 2 variables. There are 2 aspects of GLPK’s behavior for this simple problem which puzzle me:

It does not give a certificate of optimality. The termination status is “other”, and it doesn’t give an upper and lower bound as it does for other linear problems, for example, any of the problems in this workbook: https://jckantor.github.io/CBE30338/06.04-Linear-Production-Model-in-Pyomo.html

Secondly, by providing a different initial value for the decision variable, I’m able to achieve a lower objective function value than is found by GLPK. Simply change the initial value of L to be 0 instead of 1. This probably is related to the fact that there’s no certificate of optimality.

Below is a minimal example to reproduce this behavior.

Steps to reproduce the issue

def setup_model():
    m = pe.ConcreteModel()

    # add the decision variables
    def l_init(m, ix):                                                                    
        return 0                                                                      

    m.del_component('l')                                                                  
    m.add_component('l', pe.Var(range(6), within=pe.Reals, initialize=l_init))

    def L_init(m, i, i_q_k):                                                              
        return 1  # <---- change to 0 to observe lower objective function value 

    m.del_component('L')                                                                  
    m.add_component('L', pe.Var(range(3), range(6), within=pe.Reals, initialize=L_init))

    return m

m = setup_model()

m.del_component('c1')
m.add_component('c1', pe.Constraint(expr = 0 <=  m.L[0,0] + m.L[0,1] + m.L[0,2] + m.l[0]))

m.del_component('c2')
m.add_component('c2', pe.Constraint(expr = 0  <=  m.L[0,3] + m.L[0,4] + m.L[0,5]+ m.l[3]))

m.del_component('c8')
m.add_component('c8', pe.Constraint(expr = m.L[0,3] + m.L[0,4] + m.L[0,5] + m.l[3]  <=  1 + 1 - (m.l[0]  + m.L[0,0] + m.L[0,1] + m.L[0,2])))

m.del_component('c9')
m.add_component('c9', pe.Constraint(expr = m.L[1,0] + m.L[1,1] + m.L[1,2] + m.l[1]  <=  1))

# === Specify the objective function ===
m.del_component('objective_function')
m.add_component('objective_fn_expression', pe.Objective(expr = m.L[0,0] + m.l[0], sense = pe.minimize))
solver = pe.SolverFactory('glpk')
results = solver.solve(m)

Output

results
>>>> {'Problem': [{'Name': 'unknown', 'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 13, 'Number of nonzeros': 21, 'Sense': 'minimize'}], 'Solver': [{'Status': 'ok', 'Termination condition': 'other', 'Statistics': {'Branch and bound': {'Number of bounded subproblems': 0, 'Number of created subproblems': 0}}, 'Error rc': 0, 'Time': 0.009073019027709961}]}

pe.value(m.objective_fn_expression.expr)
>>>> 1

Maybe I’m missing something, but I’m not able to find anything about status "other" on the website: https://www.gnu.org/software/glpk/. Since GLPK is a linear solver, why is it having difficulty obtaining the optimal solution to a simple problem like this?

Information on your system

Pyomo version: 6.4.1 Python version: 3.10 Operating system: MacOS Monterey How Pyomo was installed (PyPI, conda, source): conda Solver (if applicable): GLPK

[makansij]

I've discovered that GLPK does not tell the user if the initial point is infeasible. Using the example above:

pe.value(m.c8.expr)
>>>> False

However, this still doesn't explain why it cannot find the optimal solution when using all 0's for the initial value of L?

I do know however, that GLPK would tell me if the problem were unbounded, since when constraints c1 and c2 are removed, it reports "unbounded".

[makansij]

It seems as if glpk also requires that there be no free variables in the optimization. In the example above m.L[1,0] + m.L[1,1] + m.L[1,2] + m.l[1] is unbounded below. Even though it doesn't play a role in the objective function, since it is unconstrained, glpk returns status other.