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

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.

[Matt J]

Please make life easy on us and attach all the Matlab variables need to run the code in a .mat file.

[Ozan Burak Ericok]

Hello Matt, I editted my question.

Answer 1

[Matt J]

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.

[Ozan Burak Ericok]

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.

[Matt J]

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.

[Ozan Burak Ericok]

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.

Answer 2

[Matt J]

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