Intlinprog can not find feasible solution even though there is one

I am currently using the 'intlinprog' solver to solve a MIP problem. In some cases, I got the following result

"No feasible folution found."

"Intlinprog stopped because no point satisfies the constraint."

However, I do have one solution, which I already testet and it satisfies the linear constraint, integer constraint, and also the lower and upper bound constraint.

Then I set this "feasible" solution as the initial point x0 for the solver. It still gave me the same message.

Afterwards, I followed the suggestion (https://support.gurobi.com/hc/en-us/community/posts/360056102052-Test-if-the-start-solution-is-feasible) and also defined the lower and upper bound of the variable exactly as this feasible solution. This means that in the feasible space there is only this point. After adding all these constraints, "intlinprog" finally found the feasible point.

Can any one please help me? Is there something that I have missed to configure??

load('Asen.mat')
load('a.mat')

% Extract the dimension of Asen, the transform matrix
% m - row number
% n - column number
m = size(Asen,1);
n = size(Asen,2);    

% Artificially generate the result 'ob'   
ob = Asen*a;

% define value of M in big M-constriant
M = 1e5;

% define the indicies of two variables in 'intlinprog'
xvar = 1:n;
yvar = n+1:2*n;
% set the intcon variable to the indicies of y -> put integer constraint on y 
intcon = yvar;
%     intcon = [];

% set lower and upper bound of both variales
% x will be positive but not greater than M, since y should be binary vector, it is constraint by bound  0<=y<=1 
lb = zeros(2*n,1);
ub = ones(2*n,1);
ub(xvar) = M;

% Construct linear equality and inequality constraints 
% The inequality constraint defined by A and b make y to be 1 when x~=0, and y to be 0 when x==0
Atemp = eye(n);
A = [Atemp, -M*Atemp];
b = zeros(n,1);
% The equality constraint defined by Aeq and beq put the constraint Asen*x==ob.
Aeq = [Asen,zeros(m,n)];
beq = ob;

% Objective function in 'inlinprog': counting the number of nonzeros(calculating the sum of y vector)
f = [zeros(n,1);ones(n,1)];

%optimize
option = optimoptions('intlinprog','IntegerTolerance',1e-6,'ConstraintTolerance',1e-3,'CutGeneration','advanced');
[xLinInt, fval, exitFlagInt, output] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,[],option);

% extract the variable x if there is a solution
if ~isempty(xLinInt)
    xval = xLinInt(xvar);
end

First a vector 'a' in length 103 is loaded. 'Asen' is a predefined matrix with dimension 202*103. Then 'ob' is obtained with 'ob=Asen*a'. The solver 'intlinprog' is supposed to find the sparsest solution 'a' satisfying 'ob=Asen*a', which is not necessarily be the original 'a' resulting 'b'. In the model, variable 'x' represent the value of 'a', binary variable 'y' in the solver represents the nonzeroness of 'x'. The objective function is to count the number of nonzero elements in x, which equals to the sum of all elements in y.

The matrix 'Asen' I used is in the attachment along with vector 'a'.

It is very strange that Matlab 2019a is not able to find a single solution. At least the original vector a I used to generate should be one feasible solution.

[Mario Malic]

People can only guess what might be the solution, if you do not post a bit of code.

Post the code where you call optimization function and its parameters at least.

[Man Sun]

Thanks. I have uploaded my code.

[Walter Roberson]

R2020a is able to handle the code, and replies with

LP:                Optimal objective value is 0.043592.                                             

Heuristics:        Found 1 solution using ZI round.                                                 
                   Upper bound is 103.000000.                                                       
                   Relative gap is 0.00%.                                                          


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the
default value). The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-06 (the selected value).

However, R2019b says that there is no solution.

[Man Sun]

Hi, I have edited my code to make it a minimal example and uploaded one more vector. Can you please run the code again in Matlab R2020a to see if it can solve the model?

[Walter Roberson]

No, it says that is not feasible.

The challenge appears to be Aeq, beq. If you zero out rows 169 onward of both, then you get a solution.

If you zap out row 106, and rows 173 onward, you get a solution. So row 106 conflicts with 169, but there is at least one more conflict involving row 173 (there might be other conflicts.)

intcon does not make a difference; A and b do not make a difference. lb and ub do not make a difference. There are internal conflicts in Aeq, beq.

[Man Sun]

It is so strange because I believe there is one feasible point x0.

x0 = [a;(a~=0)];

I checked that satisfys the constraint , , and .

Then I make the following change without changing any constraint:

• run the solver with initial solution x0

[xLinInt, fval, exitFlagInt, output] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,x0,option);

Suddenly, the model becomes feasible again. But somehow, is not able to locate this feasible point with other initial solution or empty initial solution. Could this be a problem of the solver itself?

[Walter Roberson]

I see what you mean.

I will look at this again later after I wake up.

[Walter Roberson]

Unfortunately all of the useful parts of intlinprog are inside .p files, so I cannot track to see why the difficulty is showing up.

[Man Sun]

Okay. Thanks!

Answer 1

[Matt J]

This version finds a solution. Basically, I pre-normalize max(Asen(:)) and max(a(:)) to 1 so that the constraint matrix elements don't vary so widely in order of magnitude:

load('Asen.mat')
load('a.mat')

Asen(~any(Asen,2),:)=[]; %get rid of all-zero rows

Asen=Asen/max(Asen(:)); %normalize

M0 = 1e5;
ob = Asen*a/M0; %Normalize upper bound of a(i) to 1
M=1;

m = size(Asen,1);
n = size(Asen,2);    


% define the indicies of two variables in 'intlinprog'
xvar = 1:n;
yvar = n+1:2*n;
% set the intcon variable to the indicies of y -> put integer constraint on y 
intcon = yvar;
%     intcon = [];
% set lower and upper bound of both variales
% x will be positive but not greater than M, since y should be binary vector, it is constraint by bound  0<=y<=1 
lb = zeros(2*n,1);
ub = ones(2*n,1);
ub(xvar) = M;
% Construct linear equality and inequality constraints 
% The inequality constraint defined by A and b make y to be 1 when x~=0, and y to be 0 when x==0
Atemp = eye(n);
A = [Atemp, -M*Atemp];
b = zeros(n,1);
% The equality constraint defined by Aeq and beq put the constraint Asen*x==ob.
Aeq = [Asen,zeros(m,n)];
beq = ob;
% Objective function in 'inlinprog': counting the number of nonzeros(calculating the sum of y vector)
f = [zeros(n,1);ones(n,1)];
%optimize
option = optimoptions('intlinprog');
[xLinInt, fval, exitFlagInt, output] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,[],option);



% extract the variable x if there is a solution
if ~isempty(xLinInt)
     xval = xLinInt(xvar)*M0; %Restore original scale by multiplying with M0.
end

[Man Sun]

I believe during the normalization and scaling some small values can be even smaller and will fall into the tolerance of the algorithm.

After running the above code, the xval doesn't fullfill the euqlity constraint.

>> load('Asen.mat')
load('a.mat')
>> norm(Asen*xval-Asen*a)

ans =

   14.8282

The equation Asen*x=ob actually only have one solution a because rank(Asen)=n and the number of rows m is greater than the number of columns n, basically an overdetermined but consistent system. So, the equality constraint maybe too strict for 'intlinprog' solver as there is only one point satisfying the constraint. 'intlinprog' seems more suitable for problem with large feasible space an then pick one according to the objective function. This may cause the failuer of intlinprog.

[Matt J]

The equation Asen*x=ob actually only have one solution a because rank(Asen)=n and the number of rows m is greater than the number of columns n, basically an overdetermined but consistent system...This may cause the failuer of intlinprog.

I'm sure that can happen depending on the intlinprog algorithm parameters that you set, but note that the rank is subjective. You have used the default tolerance to calculate the rank, but with a different tolerance, it will be smaller than n,

   >> n-rank(Asen,1e-6)

ans =

     6

After running the above code, the xval doesn't fullfill the euqlity constraint.

The equality constraint is fulfulled to the default tolerance,

>> norm(Aeq*xLinInt-beq)

ans =

   1.4583e-08

You can tighten the constraint tolerance if you wish, but I think it is more appropriate to assess the solution in terms of percent error. The percent error of the current solution is already quite small both according to the L2-norm and the sup-norm,

>> norm(Asen*xval-Asen*a,inf)/norm(Asen*a,inf)*100

ans =

   4.4702e-04

>> norm(Asen*xval-Asen*a)/norm(Asen*a)*100

ans =

   7.8977e-05

Answer 2

[Man Sun]

I have asked Mathworks Support. They give two nice workarounds as following:

We noticed that the constraint matrices and the problem in general will benefit from a better scaling as there continuous variables that can go upto 1000 and binary variables that can take only 0/1. So improvement on the model formulation will help "intlinprog" a lot, especially its presolve/preprocess methods.

Now, with the formulation as is, the following are two workarounds two workarounds for "intlinprog" to find the optimal.

1. Naive scaling:

[xLinInt, fval, exitFlagInt, output] = intlinprog(f,intcon,A,b,Aeq/1000,beq/1000,lb,ub);

2. Turning off LPPreprocess:

option = optimoptions('intlinprog','LPPreprocess','none');

[xLinInt, fval, exitFlagInt, output] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,option);

They both find the solution.

[Matt J]

I believe what I proposed in my 2nd answer is virtually the same as their rescaling strategy.

They both find the solution.

When I run these suggestions, I obtain

>> n-fval

ans =

    4.0000

whereas the solution you claim is the unique one should satisfy n=fval. As I mentioned here, you should not expect rank() to give you an objective prediction of what the solution space will be.

Answer 3

[Mario Malic]

I am not sure if this comment will be useful at all, since I don't understand the code completely. (Not great in optimization either, so someone else might be able to provide better comment)

For me - Line 23 throws an error.

I'll refer to the objective function, if you're trying to minimize it, what is the parameter that affects its value? The variable n doesn't change at all. So, for whatever x it is, its value will stay the same.

Suggestions that you could try:

• don't use some constraints, just replace their entry with []
• try with an initial point x0

[Man Sun]

Yes. You are correct. The value of x is of no meaning to the objective function. The objective function calculates , which calculates the sum of all binary elements in y. Since when and when , the value of the objective function is the l0 norm of x.

Answer 4

[Matt J]

Aside from the error mentioned by Mario, the line

  A = [Atemp, -M*Atemp];

is overwriting your original A matrix (the one contained in A.mat). Therefore, none of your constraints represent the condition A*a=b.

I have the vague impression that you are trying to do linearly constrained L1 norm minimization. If so, you might be interested in my minL1lin utility on the File Exchange, which already provides this functionality.

https://www.mathworks.com/matlabcentral/fileexchange/52795-minl1lin

[Man Sun]

Hi. You are correct.

I have changed the name of two variables and reran the code independently inside matlab. There should be no problem running it by now. Here I am trying to do the l0 minimization. The l1 minimization is already implemented using the 'fmincon' solver, which works really great. I really like to compare these two minimization methods.

[Matt J]

The l1 minimization is already implemented using the 'fmincon' solver, which works really great.

No, fmincon will not reliably minimize an L1-norm objective. fmincon assumes that the objective is differentiable, which the L1-norm is not.

As an example of where fmincon fails, you can try the following,

>> A=[0.9572    0.8003
    0.4854    0.1419];
>> b=A(:,1);
>> [x,fval]=fmincon(@(x) norm(A*x-b,1),[0,0].',[1 1],1)

The solution should be x=[1;0] with fval=0, but instead you get (in R2019a) the poorly converged result,

x =

    1.0112
   -0.0384


fval =

    0.0200

However, with minL1lin

>> [x,fval]=minL1lin(A,b,[1 1],1)

Optimal solution found.


x =

     1
     0


fval =

     0

[Man Sun]

Thanks! I was not aware of that! I will definitely try your solver.

In my case, all the variables are actually positive or zero. Will this make the 'fmincon' solution correct?

[Matt J]

Nope. Here's a modification of the 2D example in which we impose explicit positivity constraints. The correct solution is x=[1;1], but instead we get,

>> [x,fval]=fmincon(@(x) norm(A*x-sum(A,2),1),[0,0].',[1 1],[2],[],[],[0;0])
 

x =

    1.0478
    0.8366


fval =

    0.0850

[Man Sun]

My model is more like this:

[x,fval]=fmincon(@(x) norm(x,1),[0,0].',[],[],A,sum(A,2),[0;0])

The objective function inside is always greater or equal to zero.

[Walter Roberson]

if you are certain that the values are nonnegative then why not just sum them?

[Man Sun]

The variables are all physical value. So they can't be negative.

I thought both are expressing the same objective function. Initially, I just thought norm(x,1) will be a more simple way to write instead of

ones(1,n)*x

[Walter Roberson]

norm(, 1) is in theory not differentiable because of the abs() and it is best to avoid theoretical issues.

Have you looked at expressing using problem based optimization? sum() is supported for that.

[Man Sun]

No. I have been only using norm( ,1) inside fmincon. I will look at all the options that have been mentioned above.