MATLAB Answers

2

linprog says 'no feasible solution found' even though there is one

Asked by Ozan Burak Ericok on 2 Oct 2019
Latest activity Edited by Matt J
on 4 Oct 2019
Hello,
I am trying to solve a linear minimization problem using linprog. Variables are in the mat file. Here is the code (Briefly, minimize the sum of unknowns (weights), no inequality constraints, no upper bound, only equality constraints and lower bound):
options = optimoptions('linprog','Display','iter','Algorithm','dual-simplex',...
'Preprocess','none' );
[ weights, fval, exitflag ] = linprog( ones(num_edges,1), '', '', xAeq, xbeq,...
lb, [], options );
I am using MATLAB R2018a. When I run the code with the given variables the code says "No feasible solution found". However, I know there is a solution. For example, 'weights' in the 'data.mat' is a solution.
I changed the algorithm that I use, but none of them worked. I couldn't find the cause of this problem. I would be appriciated if you could help me with this or could suggest a different method/algorithm.
Thank you in advance.

  2 Comments

Please make life easy on us and attach all the Matlab variables need to run the code in a .mat file.
Hello Matt, I editted my question.

Sign in to comment.

2 Answers

Matt J 님의 답변 3 Oct 2019
Matt J 님이 편집함. 3 Oct 2019
 채택된 답변

Basically, the problem is that your equality constraint matrix xAeq is full rank:
>> rank(xAeq)
ans =
8
>> cond(xAeq)
ans =
10.206330042860493
I think this was unintentional on your part, because it means the constraints can have no more than one unique solution. Normally, you wouldn't make your constraints this tight, since it trivializes the optimization.
Anyway, the point is that this is apparently a numerically sensitive situation for linprog. It cannot tell the difference between the case where there is a single solution (up to floating point errors) and the case where there are none at all. It is worth pointing out though that the weights vector that you've provided in your data.mat file, and which you believe to be feasible, doesn't satisfy the constraints especially well. This is easiest to see in long precision.
>> format long; [xAeq*weights,xbeq]
ans =
0.000003546232588 0
-0.000003546232595 0
0.000000000000000 0
-0.000003546232584 0
0.000003546232600 0
-0.000003482138166 0
0.000000000000000 0
0.000003482138164 0
1.000000000000000 1.000000000000000
The least squares solution to your equality constraints also don't even satisfy the equations all that well,
>> weightsLSQ=xAeq\xbeq;
>> norm(xAeq*weightsLSQ-xbeq)
ans =
4.970015269326972e-06
I still do find it a little curious that your provided weights are considered infeasible by linprog since they do satisfy the constraints to within the ConstraintTolerance parameter. Hopefully, someone from the MathWorks will comment on that.

  3 Comments

That's exactly the problem. I found the weights by hand so there might be some numerical precision issues, but they are within the ConstraintTolerance and I think linprog should have considered it feasible. This is exactly why I posted the question.
If you want to force linprog to treat your constraints in this loose sense, one way is to replace your equality constraints with slightly relaxed inequality constraints:
load data
options = optimoptions('linprog','Display','iter','Algorithm','dual-simplex',...
'Preprocess','none');
delta=options.ConstraintTolerance/2;
options.ConstraintTolerance=delta;
A=[xAeq;-xAeq];
b=[xbeq+delta;-(xbeq-delta)];
[ weightsOptimal, fval, exitflag ] = linprog( ones(num_edges,1), A,b,[],[],...
lb, [], options );
This does find a solution, but as expected, the optimal weights are not not very different from the hypothetical weights:
Iter Time Fval Primal Infeas Dual Infeas
0 0 0.000000e+00 9.504933e-01 0.000000e+00
8 0 7.362136e+00 0.000000e+00 0.000000e+00
Optimal solution found.
>> [weights,weightsOptimal]
ans =
1.000000000000000 0.999950000000000
1.000000000000000 0.999884805203178
1.000000000000000 0.999735906204328
1.000000000000000 0.999659481819993
0.823906258837848 0.823811349123758
0.823906258837848 0.823711349123740
1.427047500981146 1.426746509993818
0.288677912917816 0.288636790167471
This had to be the case because, again, since xAeq is full column rank, the feasible set is approximately a singleton.
Hello Matt,
Sorry it took a little bit for me to test all of my cases. I was able to solve the 42 cases that I have with your approach.
Relaxing the equaility constraints slightly is essentially equivalent to what I did for practical purposes (I am okay with these approximate weights). Thank you for your help.
If I find another way of handling the equalities in the future, I will add a new comment.
Best.

Sign in to comment.


Matt J 님의 답변 4 Oct 2019
Matt J 님이 편집함. 4 Oct 2019

Here's another way you could solve it:
if rank(xAeq)==num_edges
weightsOptimal=xAeq\xbeq;
if ~all(weightsOptimal>=lb-options.ConstraintTolerance)
warning 'Solution doesn''t meet constraints within ConstraintTolerance'
end
else
weightsOptimal = linprog( ones(num_edges,1), [],[],xAeq,xbeq,...
lb, [], options );
end

  0 Comments

Sign in to comment.