When does intlinprog give wrong solutions?

Hello everyone, I am using intlinprog to solve integer programs and it seems like I am getting non-optimal solutions and I'm wondering why that might be. Specifically, for the problem

f = [1;1;1;1]
Aeq=[32,45,-41,22;-82,-13,33,-33;23,-21,12,-12]
beq=[214;712;331]
ub=[Inf;Inf;Inf;Inf]
lb=zeros(4,1)
intcon=[1:4]
x=intlinprog(f,intcon,[],[],Aeq,beq,lb,ub)

Matlab gives as the solution

x =

   1.0e+04 *

    1.3919
    4.5174
    6.9744
    1.7340

But the solution

    x=
3716
12057
18587
4582

is smaller and meets all the constraints.

Thanks for your help, Max.

[Matt J]

What MATLAB version are you using?

[Maximilian Nohr]

I am using R2015b on a Mac

Answer 1

[Matt J]

It looks like it could be a numerical stability issue. Small perturbations in the problem data can cause jump discontinuities in the solutions to integer programs, because of the inherent discontinuity of integer constraints. As a simpler example, consider the code below. In ideal , infinite precision math, the solution is x=[1;0] for any -0.1<a<0 and x=[0;10] for any 0<a<0.1. In other words, the solution has a jump discontinuity as a function of "a" at a=0.

a=+1e-6;
z=1+a;
m=-(1+2*a)/10;

A=[-1,+m];
b=-z;

f=[1;-a];
lb=[0;0];
ub=[10;10];

options=optimoptions(@intlinprog,'Display','none','TolInteger',1e-6);
[x,fval,exitflag,output]=intlinprog(f,[1,2],A,b,[],[],lb,ub,options);

Therefore, on a finite precision computer (which is to say any real-world computer), you will get unpredictable behavior for "a" small enough in magnitude. On my machine in R2015a, I get the wrong solution when I set a=-1e-8,

x =

    0.0000
   10.0000

[Maximilian Nohr]

Hey, thanks a lot for your quick answer!

[Matt J]

No problem. It might also be worth adding that intlinprog in R2015a (you didn't say what version you were using) does apparently have some bugs, and I can't rule out that they may be playing a role as well.

As an example, the modification of my toy problem below gives me non-integer results, even though intcon=1:2.

a=+.01;
z=1+a;
m=-(1+2*a)/10;

A=[-1,+m];
b=-z;

f=[1;0];
lb=[0;0];
ub=[];

options=optimoptions(@intlinprog,'Display','none');
[x,fval,exitflag,output]=intlinprog(f,[1,2],...
                         A,b,[],[],lb,ub,options);

 >> exitflag, x,fval

exitflag =

       1

x =

           0
      9.9020

fval =

       0