linprog error message: the dual appears to be infeasible (and the primal unbounded).

Hi I am trying to solve a linear optimization problem using linprog but I get the following errors most of the times (exitflag= -3 or -2).

*Exiting: One or more of the residuals, duality gap, or total relative error has stalled

*Exiting: One or more of the residuals, duality gap, or total relative error has grown 100000 times greater than its minimum value so far

I think the the problem is that some of my constraints are very small numbers compared to others. For example, A and b matrices look something like the following.

A = [1 1 0 0; 0 0 .5 1; 0 1e-8 0 1e-9]
b = [1 0 1e-12]

Moreover, when I get this error, some variables in the output of the optimization problem are negative, while I have a constraint that the variables must be non-negative.

Any idea to solve this issue would be appreciated.

Thanks

[Mehdi]

the elements of A and b I wrote are not for a real example. I just wanted to show the scales of the numbers.

The followings are an example that I got exitflag = -2.

A =[ 
0.33  -1  0  0
0.5  -1  0  0
1  -1  0  0
-1  0  0  0
1  0  0  0
0  0  1  -1
0  0  -1  0
0  0  1  0
-5.22e-07  0  -3.69e-07  0 ]
b' = [ 0 2.46e-07  1.76e-05  0  0.1  0  0  0.1 -1e-11]
c' = [0 1 0 1]

Answer 1

[Matt J]

It might be a good idea to give us the complete example so that we can run it ourselves. One thing I notice, however, is that your inequality constraints are degenerate. In particular, the second inequality

    0.5*x(3)+x(4)<=0

together with the the non-negativity constraint implies that the only feasible x are those in the subspace where x(3)=x(4)=0. However, I'm pretty sure linprog requires that the inequality constraints specify a region with positive volume in R^4. Otherwise, the problem is unstable. Small perturbations of the data can render the problem infeasible, as for example when you replace the second constraint with

     0.5*x(3)+x(4)<= -1e-20

If you know in advance that x(3)=x(4)=0 in advance, just remove them from the problem.

[Mehdi]

Thank you Matt for reply. Actually the elements of A and b I wrote in my first post are not for a real example. I just wanted to show the scales of the numbers.

Here is an example which I got exitflag = -2.

A =[ 
0.33  -1  0  0
0.5  -1  0  0
1  -1  0  0
-1  0  0  0
1  0  0  0
0  0  1  -1
0  0  -1  0
0  0  1  0
-5.22e-07  0  -3.69e-07  0
]
b' = [ 0 2.46e-07  1.76e-05  0  0.1  0  0  0.1 -1e-11]
c' = [0 1 0 1]

[Matt J]

In the examples, posted so far, there is no need for your A,b data to contain such small numbers. You could scale the final rows of A and b by 1e7 (for example), which would make the data in that row much more comparable in magnitude with the rest of the data, and without changing the meaning of the constraint.

I would also point out that it is inefficient to have rows in your A matrix with only one non-zero value, e.g., as you have in your 5th, 7th, and 8th row. Those rows are applying simple upper and lower bounds. You should be using LINPROG's lb and ub input arguments to apply such bounds.

[Matt J]

By multiplying 'b' by a big nubmer...The exit flag is not 1 but the result is correct.

When I implement my recommendations above, I get an exitflag of 1,

    A =[ 
            0.33  -1  0  0
            0.5  -1  0  0
            1  -1  0  0
            0  0  1  -1
            -5.22e-07  0  -3.69e-07  0];

    b = [ 0 2.46e-07  1.76e-05   0 -1e-11]';
    c = [0 1 0 1]';

    lb=[0, -inf,0,-inf];
    ub=[0.1,inf,0.1,inf];

    A(end,:)=A(end,:)*1e7; b(end)=b(end)*1e7;

    options = optimoptions('linprog','Algorithm','simplex','TolFun',1e-12);
    [output, value, exitflag]=linprog(c,A,b,[],[],lb,ub,[],options)

        Optimization terminated.
        output =
           1.0e-04 *
            0.1916
            0.0933
                 0
                 0
        value =
           9.3325e-06
        exitflag =
             1

[Mehdi]

Thanks a lot Matt. This solves the exitflag issue for most of the cases but for some few cases the exitflag is still -2 and the result is not correct (the output for the second variable as well as the final value is 0).

This case, for instance:

A = [     0.5           -1
            1           -1
   -0.0010438            0];
b = [0  1.3413e-05  -1e-11]';
c = [0 1]';
lb=[0, 0];
ub=[0.1, 0.1];

A(end,:)=A(end,:)*1e7; b(end)=b(end)*1e7;

options = optimoptions('linprog','Algorithm','simplex','TolFun',1e-12);
[output, value, exitflag]=linprog(c,A,b,[],[],lb,ub,[],options)

output =
   9.5807e-09
            0
value =
     0
exitflag =
    -2

I solved this problem with a different tool and I write it here which might be helpful for thise who have similar problem.

I used the CVX in MATLAB and set the cvx_precision to high. Although it is slower than linprog, it finds the exact result for this problem.

for the same input as above:

n = 2;
cvx_begin quiet 
    cvx_precision high
    variable x(n)
         minimize(c'*x)
     subject to
         A*x <=  b;
         0 <= x <= .1;
cvx_end
x 
cvx_optval
x =
 9.5807e-09
 4.7903e-09
cvx_optval =
 4.7903e-09

Answer 2

[Alan Weiss]

I cannot reproduce your result. When I entered the following, I got an answer:

A =[ 
0.5  -1  0  0  0  0
1  -1  0  0  0  0
-1  0  0  0  0  0
1  0  0  0  0  0
0  0  1  -1  0  0
0  0  -1  0  0  0
0  0  1  0  0  0
0  0  0  0  1  -1
0  0  0  0  -1  0
0  0  0  0  1  0
-1.04e-07  0  -7.03e-08  0  -8.17e-08  0
];
b = [ 0  7.94e-06  0   0.1  0  0  0.1  0  0  0.1 -1e-11]';
x = linprog(zeros(6,1),A,b)

 Optimization terminated.

 x =

    0.0945
   86.2581
    0.0856
   75.3455
    0.0856
   75.3455

This solution indeed satisfies the constraints:

A*x-b

ans =

    -86.2108
    -86.1636
     -0.0945
     -0.0055
    -75.2599
     -0.0856
     -0.0144
    -75.2599
     -0.0856
     -0.0144
     -0.0000

Alan Weiss

MATLAB mathematical toolbox documentation

[Mehdi]

Try this out with the same A abd b:

linprog([0 1 0 1 0 1]',A,b)

[Alan Weiss]

Thank you for including the reproduction steps. Indeed, the default interior-point algoririthm has trouble with that case.

Try another algorithm.

 options = optimoptions('linprog','Algorithm','dual-simplex');
 x = linprog(f,A,b,[],[],[],[],[],options)

 Optimal solution found.

 x =

     0
     0
     0
     0
     0
     0

 A*x-b

 ans =

         0
   -0.0000
         0
   -0.1000
         0
         0
   -0.1000
         0
         0
   -0.1000
    0.0000

Or, if you don't have a recent enough version of the toolbox to have the 'dual-simplex' algorithm, try setting the Algorithm option to 'simplex' to get the same result.

Alan Weiss

MATLAB mathematical toolbox documentation

[Mehdi]

Thanks for your help. I do not have the latest version of the toolbox so I use the 'simplex' instead of 'dual-simplex'.

with simplex, I get correct result sometimes by setting ’TolFun’ to a very small number. For some cases, I still cannot get the correct result. I know that the correct result is in the order of'1e-5 to 1e-10

The case that I can get a correct result:

A =[ 
0.5  -1  0  0  0  0
1  -1  0  0  0  0
-1  0  0  0  0  0
1  0  0  0  0  0
0  0  1  -1  0  0
0  0  -1  0  0  0
0  0  1  0  0  0
0  0  0  0  1  -1
0  0  0  0  -1  0
0  0  0  0  1  0
-1.04e-07  0  -7.03e-08  0  -8.17e-08  0];
b = [ 0  7.94e-06  0   0.1  0  0  0.1  0  0  0.1 -1e-11]';
c = [0 1 0 1 0 1]';
options = optimoptions('linprog','Algorithm','simplex','TolFun',1e-12);
[output, value, exitflag]=linprog(c,A,b,[],[],[],[],[],options)
Optimization terminated.
output =
   9.6154e-05
   8.8214e-05
            0
            0
            0
            0
value =
   8.8214e-05
exitflag =
     1

But same settings do not work with the following case (notice that the exitflag is 1):

A =[ 
0.33  -1  0  0
0.5  -1  0  0
1  -1  0  0
-1  0  0  0
1  0  0  0
0  0  1  -1
0  0  -1  0
0  0  1  0
-5.22e-07  0  -3.69e-07  0]
b = [ 0 2.46e-07  1.76e-05  0  0.1  0  0  0.1 -1e-11]';
c = [0 1 0 1]';
options = optimoptions('linprog','Algorithm','simplex','TolFun',1e-12);
[output, value, exitflag]=linprog(c,A,b,[],[],[],[],[],options)
Optimization terminated.
output =
  -2.1176e-22
   3.3881e-21
            0
            0
value =
   3.3881e-21
exitflag =
     1

The above result is not correct because I have the non-negative constraints for the values.

If I do not set any values for 'Tolfun', the linprog ignores the small numbers which is not suitable for me.

options = optimoptions('linprog','Algorithm','simplex');
[output, value, exitflag]=linprog(c,A,b,[],[],[],[],[],options)
Optimization terminated.
output =
     0
     0
     0
     0
     0
     0
value =
     0
exitflag =
     1

[Matt J]

The above result is not correct because I have the non-negative constraints for the values.

The constraints are still satisfied within the TolCon parameter.

[Mehdi]

The constraints are still satisfied within the TolCon parameter.

How can I see this? As far as I understood TolCon is for GA. How can I adjust this parameter for linprog?

In fact based on my problem the result for 'value' must not become less than 1e-11 but in the above examples, it becomes 0 and 3.39e-21 for different settings which are not correct answers for my problem.

[Matt J]

Hmmm, it looks like only the dual-simplex algorithm has a TolCon paramater, see LINPROG Options.

Anyway, the point is to remember that you're working on finite precision computers and so exact math cannot be expected of them. You must expect that the constraints may be violated by small amounts due to floating point errors.

In a linear program, in particular, the solution is always at the boundary of the feasible region, so floating point errors can push the solver slightly outside of this boundary.

[Mehdi]

OK. So it seems that I have to use dual-simplex !

Do you know any tricks to get rid of that precision issue and solve the problem by simplex?

Actually I use a tricky way and apparently it solves the issue. I multiply 'b' by a big number like '1e10' and solve the problem and at the end divide the result by the same number I multiplied at the beginning. It works but still the exitflag is not 1.

See this:

A =[ 
0.33  -1  0  0
0.5  -1  0  0
1  -1  0  0
-1  0  0  0
1  0  0  0
0  0  1  -1
0  0  -1  0
0  0  1  0
-5.22e-07  0  -3.69e-07  0];
b = [ 0 2.46e-07  1.76e-05  0  0.1  0  0  0.1 -1e-11]';
c = [0 1 0 1]';
options = optimoptions('linprog','Algorithm','simplex','TolFun',1e-12);
[output, value, exitflag]=linprog(c,A,b,[],[],[],[],[],options)
Optimization terminated.
output =
  -2.1176e-22
   3.3881e-21
            0
            0
value =
   3.3881e-21
exitflag =
     1

Notice that the exitflag= 1 but the result is not correct, at least for my problem.

By multiplying 'b' by a big nubmer:

output, value, exitflag]=linprog(c,A,1e10*b,[],[],[],[],[],options);
Exiting: The constraints are overly stringent; no feasible starting point found.
output/1e10
ans =
   1.9157e-05
   9.3325e-06
            0
            0
exitflag =
    -2

The exit flag is not 1 but the result is correct.

I am not sure if this is a reliable solution! :|

[Torsten]

If you multiply b by a factor of 1e10, you will also have to multiply A with this big number to get an equivalent problem to the one you started with.

Best wishes

Torsten.

[Matt J]

In addition to what Torsten said, scaling the entire A,b data will not fix the original problem because the relative magnitudes of the last row to the other rows remains the same.

That is why I proposed scaling the final rows only, as in my Comment here,

http://www.mathworks.com/matlabcentral/answers/225887#comment_295767