Solver hangs on a trivially unbounded LP

[LucasBoTang]

Problem:

minimize x1 + x2
subject to x1 + 2 x2 = 5
3 x1 + 2 x2 <= 8
(no variable bounds)

Summary

This very small LP is mathematically primal-unbounded, but the solver does not terminate (keeps running for a long time).

Reproducible example

C interface

#include "cupdlpx/interface.h"
#include <stdio.h>
#include <math.h>

int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, (no variable bounds)
     //
     // minimize   x1 + x2
     // subject to 1 x1 + 2 x2 = 5
     //            3 x1 + 2 x2 <= 8
     // (no variable bounds)

    int m = 2; // number of constraints
    int n = 2; // number of variables

    // objective c
    double c[2] = {1.0, 1.0};

    // A as a dense matrix
    double A[2][2] = {
        {1.0, 2.0},
        {3.0, 2.0}
    };

    // describe A as dense
    matrix_desc_t A_dense;
    A_dense.m = m;
    A_dense.n = n;
    A_dense.fmt = matrix_dense;
    A_dense.zero_tolerance = 0.0;
    A_dense.data.dense.A = &A[0][0];

    // constraint bounds l <= A x <= u
    double l[2] = {5.0, -INFINITY};
    double u[2] = {5.0,  8.0     };

    // variable bounds
    double *lb = malloc(2 * sizeof *lb);
    lb[0] = -INFINITY;
    lb[1] = -INFINITY;
    double* ub = NULL;

    // Build problem
    lp_problem_t* prob = create_lp_problem(
        &A_dense,     // constraint matrix A
        c,            // objective_c
        NULL,         // objective constant
        lb,           // variable lower bounds
        ub,           // variable upper bounds
        l,            // constraint lower bounds
        u             // constraint upper bounds
    );

    // solve
    cupdlpx_result_t* res = solve_lp_problem(prob, NULL);
    lp_problem_free(prob);
}

Python Interface

import numpy as np
from cupdlpx import Model

# Problem:
#   minimize   x1 + x2
#   subject to x1 + 2 x2 = 5
#              3 x1 + 2 x2 <= 8
#   (no variable bounds)

# setup model
c = np.array([1.0, 1.0])
A = np.array([[1.0, 2.0],
              [3.0, 2.0],
            ])
l = np.array([5.0, -np.inf])
u = np.array([5.0, 8.0])
lb = np.array([-np.inf, -np.inf])
ub = None
m = Model(c, A, l, u, lb, ub)

# solve
m.optimize()

Additional note

I also tried a higher feasibility tolerance; the solver still does not terminate on this unbounded instance.

m.set_param("eps_feasibility", 1e-6) 

[ZedongPeng]

Closed by #56