Linprogr in Matlab not finding a solution when a solution exists

I am solving a linear programming in Matlab using linprogr. All the matrices are attached.

clear
rng default
load matrices
f=zeros(size(x,1),1);
possible_solution = linprog(f,Aineq,bineq,Aeq,beq, lb, ub); 

The output is

No feasible solution found.
Linprog stopped because no point satisfies the constraints.

From some theory behind the problem, I know that the vector x (uploaded together with the other matrices) should be a solution of linprog(f,Aineq,bineq,Aeq,beq, lb, ub).

In fact, x seems to satisfy all the constraints, except the equality constraints. However, regarding the equality constraints, the maximum absolute difference between Aeq*x and beq is around 3*10^(-15). Hence, also the equality constraints are almost satisfied by x.

load x
check_lb=sum(x>=0)==size(x,1);
check_ub=sum(x<=1)==size(x,1);
check_ineq=(sum(Aineq*x<=bineq)==size(Aineq,1));
check_eq=max(abs((Aeq*x-beq)));

Therefore, my question is: why linprog(f,Aineq,bineq,Aeq,beq, lb, ub) does not find x as solution? It does not seem to be a problem of the equality constraints, because if I remove the inequality constraints and run

rng default
possible_solution_no_inequalities = linprog(f,[],[],Aeq,beq, lb, ub); 

the program finds a solution. Is it a matter of numerical precision? How can I control for that?

[Bruno Luong]

Four year (sept 2023) later where the problem is reported and still the same issue with linprog.

[Bruno Luong]

UPDATE R2024A has a new default algorithm 'dual-simplex-highs' and this issue seems to be somewhat resolve

load matrices
f=zeros(size(x,1),1);
possible_solution = linprog(f,Aineq,bineq,Aeq,beq, lb, ub); 

Optimal solution found.

Sill fail on @Simão Pedro da Graça Oliveira Marto example below.

Answer 1

[Bruno Luong]

It looks to me LINPROG is flawed or buggy in your case. I try to remove suspected constraints and change beq so that the equality is satisfied by x, and LINPROG still fails.

load matrices
f=zeros(size(x,1),1);
Aeq = round(Aeq);
beq = Aeq*x;
keep = max(abs(Aineq),[],2)>1e-10;

all(x>=lb)

ans = logical
   1

all(x<=ub)

ans = logical
   1

all(Aineq*x<=bineq)

ans = logical
   1

all(Aeq*x==beq)

ans = logical
   1

% This will fail to return solution
options = optimoptions('linprog','Algorithm','interior-point','ConstraintTolerance',1e-3);
xlinprog = linprog(f,Aineq(keep,:),bineq(keep,:),Aeq,beq, lb, ub, options); 

The problem is infeasible.

Linprog stopped because no point satisfies the constraints.

But if I call QUADPROG it will returns a solution

% This will returns a solution
H = sparse(size(x,1),size(x,1));
xquadprog = quadprog(H,f,Aineq(keep,:),bineq(keep,:),Aeq,beq, lb, ub);

Warning: Quadratic matrix (H) is empty or all zero. Your problem is linear. Consider calling LINPROG instead.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

[CT]

Thanks. I tried to call gurobi from Matlab (in place of linprogr): gurobi always finds a solution.

[Rub Ron]

I am very much dissapointed and upset by linprog and Matlab. I had spend a lot of time dealling with errors in my algorithm. I was putting the blame in some flaw in my algorithm and precomputation of variables, but at the end the issue was with linprog. I switched now to gurobi and my algorithm does perform as expected. It was a huge waste of time and effort, and the documentation on linprog is poor, it does not even flag cases with "challlenging" data.

Answer 2

[Derya]

Hello CT,

linprog with 'Preprocess' option set to 'none' will give you a solution. Then, by decreasing the 'ConstraintTolerance' you may get a solution with better feasibility measures. Since the problem is numerically challenging(*) you may need to examine any solution found (by any solver) carefully.

Thank you very much for providing the example. Using it, we'll try to improve linprog so that it can solve these type of problems with its default settings.

Kind Regards,

Derya

(*)

1. there are very small coefficient in matrices which are below some of the tolerances used in our numerical algorithms.

2. the ratio between the largest and smallest absolute values of coefficients in the constraint matrices is large.

[Simão Pedro da Graça Oliveira Marto]

Hi Derya,

I found a simpler problem that has the same issue. In only has 99 variables and 1 linear constraint. The relevant arrays are attached. The following code runs the problem:

load('linprog_failure_case', 'obj', 'Aeq', 'beqMin', 'beqMax', 'lb', 'ub');

options = optimoptions("linprog", "Display", "iter");
[wMin, A0] = linprog(+obj,[],[],Aeq,beqMin,lb,ub,options);
[wMax, B0] = linprog(-obj,[],[],Aeq,beqMax,lb,ub,options);

x = fzero(@(x)Aeq*(lb+x*(ub-lb))-beqMin, 0.5);
w = lb+x*(ub-lb); % This point solves constraints for first problem

assert(all(w>=lb) && all(w<=ub) && abs(Aeq*w-beqMin)<1e-15)

Both linprog return the same error: "Linprog stopped because no point satisfies the constraints."

However, the variable w meets all constraints for the first problem, and a similar solution can be obtained for the other one.

Setting the option "Preprocess" to "none" also solves this problem, I am only commenting to give you another example, in case it is useful for you.

The difference between upper and lower bounds is very small, between the order of 1e-12 and 1e-8. But when I randomly generate problems with the same feature, linprog seems to always work, so there may be more to it.

[Derya]

Thank you very much for providing the problem instance, Simão Pedro.

Kind Regards,

Derya

[Bruno Luong]

@Simão Pedro da Graça Oliveira Marto and @Derya

The issue of problem submited by Semao seems to be different than CT. If we set

options = optimoptions('linprog','Algorithm','interior-point')

In Simao problem linprog will find solution. This is not the case with CT's problem.

quadprog works for both problems.

Answer 3

[Matt J]

The problem seems to be that the inequality constrained region has barely any intersection with the equality constrained region. This means that the feasible set, if it exists, is very small, making the solution very unstable numerically, as well as difficult for linprog to find.

As evidence of this, the code below shows that by enlarging the inequality constrained region very slightly to Aineq*x<=bineq+1e-10, a solution is found. Note that this is after a pre-normalization of the rows of the constraint matrices.

load('matrices.mat')

[Aineq,bineq]=normalizeConstraints(Aineq,bineq,'<=');
[Aeq,beq]=normalizeConstraints(Aeq,beq,'==');

opts=optimoptions('linprog','ConstraintTol',1e-9);
f=zeros(size(x,1),1);
xp = linprog(f,Aineq,bineq+1e-10,Aeq,beq, lb, ub,opts); 

Optimal solution found.

function [A,b]=normalizeConstraints(A,b,op)

    idx=~any(A,2); %find rows with all zeros

    switch op

        case '=='

        assert( all(b(idx)==0) , 'Trivially infeasible')

        case '<='

        assert( all(b(idx)>=0) , 'Trivially infeasible')

    end

    A(idx,:)=[]; %discard rows with all zeros
    b(idx)=[];
    
    s=vecnorm(A,2,2);
    
    A=A./s; %normalize rows of A to unit l2-norm 
    b=b./s;

end

[Bruno Luong]

As in my post quadprog ables to find a feasible solution, with the same constraints.

IMO there is no reason why different optimization functions not having to deal with non-linear constraints (namely linprog, quadprog, lsqlin, ...) cannot share the same subfunction that search of feasible solution at the preprocessing step.

[Matt J]

As in my post quadprog ables to find a feasible solution, with the same constraints.

Yes, but what meaning does that solution have if the feasibility of the problem is fundamentally unstable, as my post conjectures?

IMO there is no reason why different optimization functions not having to deal with non-linear constraints cannot share the same subfunction that search of feasible solution at the preprocessing step.

They could, but is it the best thing? The default linear programming algorithm works in dual space, whereas quadprog does not. Would it be better for linprog to work in primal space in order to be unified with quadprog? In a linear program, you can get strong information from the dual. If the dual problem appears unbounded you can conclude that the problem is infeasible and not search for an initial feasible point at all.

[Bruno Luong]

The last time I check I has an impression that the preprocessing step for both linprog and quadprog use some heuristic method based on linprog to find the feasible solution.I don't care what is behind the scene, whereas it uses dual or primal formuation (actually I do care). All I need is it is robust and reliable.

I just point out that TMW can use both quadprog and linprog to analyze the issue here and make a more robust preprocessing for all functions.

[Bruno Luong]

@Matt J "by enlarging the inequality constrained region very slightly to Aineq*x<=bineq+1e-10"

After the normalization by row norm of Aineq a big junk of the bineq is still as low as 1e-12. So actually the addition 1e-10 value is not slighly IMO. The expected unit of the solution x i snot known and it should be somehow play a role too.

[Matt J]

All I need is it is robust and reliable.

Well, it can't be if the feasible set is unstable. Imagine the feasible set was just one point. Obviously small perturbations in the constraint matrix data can render the feasible set empty, and no initialization algorithm can give meaningful results.

The expected unit of the solution x i snot known and it should be somehow play a role too.

An increment of 1e-10 in bineq means that the L2-distance to all the constraint boundaries has increased by 1e-10, a percent change in the given of 4.2e-9%.

You could still argue that the units of the different x(i) may be poorly scaled relative to one another, but if so that wouldn't be linprog's fault. linprog has to assume that the relatve scaling is well behaved.

[Bruno Luong]

I'm not convince by your normalization. You seems to focus all on 2-norm of rows of Aineq,and normlize wrt it.

What I see is the original bineq has constant value of 1e-10. That indicates me that the units of bineq elements are coherent and one should not touch it by scaling it non uniformly.

When you perform normalization and add 1e-10 (after normalization), you actually multiply the rhs of half of the inequalities by 66 to 214 from original value of bineq. This is not a small modification of the problem by all mean.

Might be quadprog and gurobi succeed to fid feasible solutions and linprog fails for this specific example is just pure (un)luck.

But sorry Matt, but if I'm the person who is responsible for linprog implementation, I would NOT blame user setup with such quick analysis as you did.

[Matt J]

You seems to focus all on 2-norm of rows of Aineq,and normlize wrt it. That indicates me that the units of bineq elements are coherent and one should not touch it by scaling it non uniformly.

The normalization doesn't change the region identified by the constraints. It only makes it so that bineq(i) measures the distance of the constraint boundary Aineq(i,:) from the origin in the same units as . If you don't normalize, there is no way to assess the positions of the boundaries in R^n w.r.t. each other.

Might be quadprog and gurobi succeed to fid feasible solutions and linprog fails for this specific example is just pure (un)luck.

Yes!

But sorry Matt, but if I'm the person who is responsible for linprog implementation, I would NOT blame user setup with such quick analysis as you did.

You have to put some responsibility on the user to choose judicious units for the x(i). Suppose, due to a crazy choice of units, the user's constraints were just box bounds in R^2 with length 1 and width 1e-31. This supplies linprog with a constraint set with a volume in R^2 that is numerically indistinguishable from 0. Should linprog be expected to handle that?

[Bruno Luong]

If you don't normalize, there is no way to assess the positions of the boundaries in R^n w.r.t. each other.

I'm not protesting about the normalization by itself what I protest is that you you add 1e-10 uniformly after normalization, and claim that is small modification. This does not convince me.

You have to put some responsibility on the user to choose judicious units.

On my view the units of the example is fine, bineq has uniform value so same unit. The uni of x is fine (as you told earlier), With this assumptions the unit of elements of Aineq must be correct (uniform).

We don't know the physics and the unit of Aineq. You seem to imply that OP doesn't setup correctly. Again you focus on row norm of Aineq. I tell that it is just an assumption that might and might not be correct.

Yes!

OP wrote: gurobi always finds a solution.

gurobi has very very good luck then.

[Matt J]

Alright. How about this for additional evidence. In the modification below, all I have done is translate the feasible set, so that the given feasible x moves to 0. Thus x=0 is now a feasible point. The solver is able to determine that the problem is feasible in this case. However, no matter how I randomly vary the objective f, the solver cannot find a different optimal point xnew. This has to be because the solver sees the feasible set as a singleton, wherein x is the only feasible point.

load('matrices.mat')

[Aineq,bineq]=normalizeConstraints(Aineq,bineq,'<=');
[Aeq,beq]=normalizeConstraints(Aeq,beq,'==');


opts=optimoptions('linprog','Algorithm','interior-point');


deltaX=nan(1,10);

for i=1:numel(deltaX)
    
f=randn(size(x,1),1);

 xnew = linprog(f,Aineq,bineq-Aineq*x,Aeq,beq-Aeq*x, lb, ub,opts)+x; 
 deltaX(i)=max(abs(xnew-x));
end

Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.


Solution found during presolve.

maxChange = max(deltaX)

maxChange = 0

function [A,b]=normalizeConstraints(A,b,op)

    idx=~any(A,2); %find rows with all zeros

    switch op

        case '=='

        assert( all(b(idx)==0) , 'Trivially infeasible')

        case '<='

        assert( all(b(idx)>=0) , 'Trivially infeasible')

    end

    A(idx,:)=[]; %discard rows with all zeros
    b(idx)=[];
    
    s=vecnorm(A,2,2);
    
    A=A./s; %normalize rows of A to unit l2-norm 
    b=b./s;

end

[Bruno Luong]

Interesting. I have to think about it more tomorrow. Now it is too late and I'm lazy to dig into your code.

[Bruno Luong]

@Matt J The last code with shifted solution by x wrong since you forgot to change accordingly lb and ub. IMO the correct command is

xnew = linprog(f,Aineq,bineq-Aineq*x,Aeq,beq-Aeq*x, lb-x, ub-x,opts)+x;