Linprogr in Matlab not finding a solution when a solution exists

Question

CT 2019-8-14

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476153-linprogr-in-matlab-not-finding-a-solution-when-a-solution-exists

编辑： Bruno Luong 2024-3-26

采纳的回答： Bruno Luong

matrices.mat

在 MATLAB Online 中打开

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?

2 个评论
显示无隐藏无

Bruno Luong 2023-9-25

编辑：Bruno Luong 2023-9-25

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

Bruno Luong 2024-3-26

编辑：Bruno Luong 2024-3-26

在 MATLAB Online 中打开

matrices.mat

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 2019-8-14

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476153-linprogr-in-matlab-not-finding-a-solution-when-a-solution-exists#answer_387613

编辑：Bruno Luong 2019-8-14

在 MATLAB Online 中打开

matrices.mat

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.

2 个评论
显示无隐藏无

CT 2019-8-14

编辑：CT 2019-8-14

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

Rub Ron 2023-10-24

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 2019-8-15

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476153-linprogr-in-matlab-not-finding-a-solution-when-a-solution-exists#answer_387759

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Simão Pedro da Graça Oliveira Marto 2023-9-25

在 MATLAB Online 中打开

linprog_failure_case.mat

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 2023-9-25

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

Kind Regards,

Derya

Bruno Luong 2023-9-27

编辑：Bruno Luong 2023-9-27

在 MATLAB Online 中打开

@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 2023-9-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476153-linprogr-in-matlab-not-finding-a-solution-when-a-solution-exists#answer_1318072

编辑：Matt J 2023-9-26

在 MATLAB Online 中打开

matrices.mat

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

11 个评论
显示 9更早的评论隐藏 9更早的评论

Matt J 2023-9-26

编辑：Matt J 2023-9-26

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?

Matt J 2023-9-26

编辑：Matt J 2023-9-26

在 MATLAB Online 中打开

matrices.mat

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