intlinprog.m


% Copyright 2017, Gurobi Optimization, Inc.

function [x,fval,exitflag,output] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,x0,options)
%INTLINPROG A mixed integer programming (MIP) example using the
%   Gurobi MATLAB interface
%
%   This example is based on the intlinprog interface defined in the
%   MATLAB Optimization Toolbox. The Optimization Toolbox
%   is a registered trademark of The MathWorks, Inc.
%
%   x = INTLINPROG(f,intcon,A,b) solves the MIP problem:
%
%     minimize     f'*x
%     subject to    A*x <= b,
%                     x(j) integer, where j is in the vector
%                                   intcon of integer constraints.
%
%   For large problems, you can pass A as a sparse matrix and b as a
%   sparse vector.
%
%   x = INTLINPROG(f,intcon,A,b,Aeq,beq) solves the MIP problem:
%
%     minimize     f'*x
%     subject to    A*x <= b,
%                 Aeq*x == beq,
%                     x(j) integer, where j is in the vector
%                                   intcon of integer constraints.
%
%   For large problems, you can pass Aeq as a sparse matrix and beq as a
%   sparse vector. You can set A=[] and b=[] if no inequalities exist.
%
%   x = INTLINPROG(f,intcon,A,b,Aeq,beq,lb,ub) solves the MIP problem:
%
%     minimize     f'*x
%     subject to    A*x <= b,
%                 Aeq*x == beq,
%           lb <=     x <= ub,
%                     x(j) integer, where j is in the vector
%                                   intcon of integer constraints.
%
%   You can set lb(j) = -inf, if x(j) has no lower bound, and ub(j) = inf,
%   if x(j) has no upper bound. You can set Aeq=[] and beq=[] if no
%   equalities exist.
%
%   x = INTLINPROG(f,intcon,A,b,Aeq,beq,lb,ub,X0) solves the problem above
%   with MIP start set to X0.
%
%   You can set lb=[] or ub=[] if no bounds exist.
%
%   x = INTLINPROG(f,intcon,A,b,Aeq,beq,lb,ub,x0,OPTIONS) solves the
%   problem above given the specified OPTIONS. Only a subset of possible
%   options have any effect:
%
%     OPTIONS.Display               'off' or 'none' disables output,
%     OPTIONS.MaxTime               time limit in seconds,
%     OPTIONS.MaxFeasiblePoints     MIP feasible solution limit,
%     OPTIONS.RelativeGapTolerance  relative MIP optimality gap,
%     OPTIONS.AbsoluteGapTolerance  absolute MIP optimality gap.
%
%   x = INTLINPROG(PROBLEM) solves PROBLEM, which is a structure that must
%   have solver name 'intlinprog' in PROBLEM.solver. You can also specify
%   any of the input arguments above using fields PROBLEM.f, PROBLEM.A, ...
%
%   [x,fval] = INTLINPROG(f,intcon,A,b) returns the objective value at the
%   solution. That is, fval = f'*x.
%
%   [x,fval,exitflag] = INTLINPROG(f,intcon,A,b) returns an exitflag
%   containing the status of the optimization. The values for exitflag and
%   the corresponding status codes are:
%
%      2  stopped prematurely, integer feasible point found
%      1  converged to a solution
%      0  stopped prematurely, no integer feasible point found
%     -2  no feasible point found
%     -3  problem is unbounded
%
%   [x,fval,exitflag,OUTPUT] = INTLINPROG(f,intcon,A,b) returns information
%   about the optimization. OUTPUT is a structure with the following fields:
%
%     OUTPUT.message          Gurobi status code
%     OUTPUT.relativegap      relative MIP optimality gap
%     OUTPUT.absolutegap      absolute MIP optimality gap
%     OUTPUT.numnodes         number of branch-and-cut nodes explored
%     OUTPUT.constrviolation  maximum violation for constraints and bounds
%

% Initialize missing arguments
if nargin == 1
    if isa(f,'struct') && isfield(f,'solver') && strcmpi(f.solver,'intlinprog')
        [f,intcon,A,b,Aeq,beq,lb,ub,x0,options] = probstruct2args(f);
    else
        error('PROBLEM should be a structure with valid fields');
    end
elseif nargin < 4 || nargin > 10
    error('INTLINPROG: the number of input arguments is wrong');
elseif nargin < 10
    options = struct();
    if nargin == 9
        if isa(x0,'struct') || isa(x0,'optim.options.SolverOptions')
            options = x0; % x0 was omitted and options were passed instead
            x0 = [];
        end
    else
        x0 = [];
        if nargin < 8
            ub = [];
            if nargin < 7
                lb = [];
                if nargin < 6
                    beq = [];
                    if nargin < 5
                        Aeq = [];
                    end
                end
            end
        end
    end
end

%build Gurobi model
model.obj = f;
model.A = [sparse(A); sparse(Aeq)]; % A must be sparse
n = size(model.A, 2);
model.vtype = repmat('C', n, 1);
model.vtype(intcon) = 'I';
model.sense = [repmat('<',size(A,1),1); repmat('=',size(Aeq,1),1)];
model.rhs = full([b(:); beq(:)]); % rhs must be dense
if ~isempty(x0)
    model.start = x0;
end
if ~isempty(lb)
    model.lb = lb;
else
    model.lb = -inf(n,1); % default lb for MATLAB is -inf
end
if ~isempty(ub)
    model.ub = ub;
end

% Extract relevant Gurobi parameters from (subset of) options
params = struct();

if isfield(options,'Display') || isa(options,'optim.options.SolverOptions')
    if any(strcmp(options.Display,{'off','none'}))
        params.OutputFlag = 0;
    end
end

if isfield(options,'MaxTime') || isa(options,'optim.options.SolverOptions')
    params.TimeLimit = options.MaxTime;
end

if isfield(options,'MaxFeasiblePoints') ...
        || isa(options,'optim.options.SolverOptions')
    params.SolutionLimit = options.MaxFeasiblePoints;
end

if isfield(options,'RelativeGapTolerance') ...
        || isa(options,'optim.options.SolverOptions')
    params.MIPGap = options.RelativeGapTolerance;
end

if isfield(options,'AbsoluteGapTolerance') ...
        || isa(options,'optim.options.SolverOptions')
    params.MIPGapAbs = options.AbsoluteGapTolerance;
end

% Solve model with Gurobi
result = gurobi(model, params);

% Resolve model if status is INF_OR_UNBD
if strcmp(result.status,'INF_OR_UNBD')
    params.DualReductions = 0;
    warning('Infeasible or unbounded, resolve without dual reductions to determine...');
    result = gurobi(model,params);
end

% Collect results
x = [];
output.message = result.status;
output.relativegap = [];
output.absolutegap = [];
output.numnodes = result.nodecount;
output.constrviolation = [];

if isfield(result,'x')
    x = result.x;
    if nargout > 3
        slack = model.A*x-model.rhs;
        violA = slack(1:size(A,1));
        violAeq = norm(slack((size(A,1)+1):end),inf);
        viollb = model.lb(:)-x;
        violub = 0;
        if isfield(model,'ub')
            violub = x-model.ub(:);
        end
        output.constrviolation = max([0; violA; violAeq; viollb; violub]);
    end
end

fval = [];

if isfield(result,'objval')
    fval = result.objval;
    if nargout > 3 && numel(intcon) > 0
        U = fval;
        L = result.objbound;
        output.relativegap = 100*(U-L)/(abs(U)+1);
        output.absolutegap = U-L;
    end
end

if strcmp(result.status, 'OPTIMAL')
    exitflag = 1;
elseif strcmp(result.status, 'INFEASIBLE') ...
        || strcmp(result.status, 'CUTOFF')
    exitflag = -2;
elseif strcmp(result.status, 'UNBOUNDED')
    exitflag = -3;
elseif isfield(result, 'x')
    exitflag = 2;
else
    exitflag = 0;
end

% Local Functions =========================================================

function [f,intcon,A,b,Aeq,beq,lb,ub,x0,options] = probstruct2args(s)
%PROBSTRUCT2ARGS Get problem structure fields ([] is returned when missing)

f = getstructfield(s,'f');
intcon = getstructfield(s,'intcon');
A = getstructfield(s,'Aineq');
b = getstructfield(s,'bineq');
Aeq = getstructfield(s,'Aeq');
beq = getstructfield(s,'beq');
lb = getstructfield(s,'lb');
ub = getstructfield(s,'ub');
x0 = getstructfield(s,'x0');
options = getstructfield(s,'options');

function f = getstructfield(s,field)
%GETSTRUCTFIELD Get structure field ([] is returned when missing)

if isfield(s,field)
    f = getfield(s,field);
else
    f = [];
end

Try Gurobi for Free

Choose the evaluation license that fits you best, and start working with our Expert Team for technical guidance and support.

Evaluation License
Get a free, full-featured license of the Gurobi Optimizer to experience the performance, support, benchmarking and tuning services we provide as part of our product offering.
Academic License
Gurobi supports the teaching and use of optimization within academic institutions. We offer free, full-featured copies of Gurobi for use in class, and for research.
Cloud Trial

Request free trial hours, so you can see how quickly and easily a model can be solved on the cloud.

Search