Mixed-Integer Linear Programming (MILP) Algorithms
Mixed-Integer Linear Programming Definition
A mixed-integer linear program (MILP) is a problem with
Linear objective function, fTx, where f is a column vector of constants, and x is the column vector of unknowns
Bounds and linear constraints, but no nonlinear constraints (for definitions, see Write Constraints)
Restrictions on some components of x to have integer values
In mathematical terms, given vectors f, lb,
and ub, matrices A and Aeq,
corresponding vectors b and beq, and a set of
indices intcon, find a vector x to
solve
HiGHS MILP Algorithm
Overview of HiGHS
The intlinprog algorithm is based on the HiGHS
open-source software. intlinprog converts MATLAB®-formatted inputs and options into equivalent HiGHS arguments, and
converts the returned solution into standard MATLAB format as well.
Algorithm Outline
The algorithm performs these steps.
Get a branch/bound node (if none then stop)
Repeat until a stopping condition:
Search until a stopping condition:
Perform Plunge/Diving
Propagate domain
Prune nodes, update bounds, exit if infeasibility detected or
ObjectiveCutOffoption value is reachedIf restart conditions are met, return to step 2
Install the next node:
Choose and evaluate a node
If evaluation prunes this node, return to step 5
Generate cuts for the node
If domain is infeasible, cut off the node, and bring open nodes into the node queue
Update the basis
Go to step 4
Presolve
Usually, it is possible to reduce the number of variables in the problem (the number of components of x), and reduce the number of linear constraints. While performing these reductions can take time for the solver, they usually lower the overall time to solution, and can make larger problems solvable. The algorithms can make solution more numerically stable. Furthermore, these algorithms can sometimes detect an infeasible problem.
Presolve steps aim to eliminate redundant variables and constraints, improve the scaling of the model and sparsity of the constraint matrix, strengthen the bounds on variables, and detect the primal and dual infeasibility of the model. For background, see Andersen and Andersen [3] and Mészáros and Suhl [10].
During mixed-integer program preprocessing, intlinprog
analyzes the linear inequalities A*x ≤ b along with
integrality restrictions to determine whether:
The problem is infeasible.
Some bounds can be tightened.
Some inequalities are redundant, so can be ignored or removed.
Some inequalities can be strengthened.
Some integer variables can be fixed.
For background about integer preprocessing, see Achterberg et al. [1].
Evaluate Root Node
To evaluate the root node, the algorithm performs the following steps.
Detect symmetry and simplify problem.
Evaluate root LP.
Generate and add LP cuts (see Cut Generation).
Apply randomized rounding.
Generate and add more LP cuts.
Perform cut generation and heuristics in a loop.
Check for restart conditions, restart at step 2 if warranted.
After the root node evaluation completes, the algorithm proceeds with Branch-and-Bound.
Cut Generation
Cuts are additional linear inequality constraints that
intlinprog adds to the problem. These inequalities
attempt to restrict the feasible region of the LP relaxations so that their
solutions are closer to integers. For background about cut generation algorithms
(also called cutting plane methods), see Cornuéjols [6] and Atamtürk, Nemhauser, and Savelsbergh [4].
Plunge/Diving
To find integer-feasible points, intlinprog uses
heuristics that are similar to branch-and-bound steps, but follow just one
branch of the tree down, without creating the other branches. This single branch
leads to a fast “dive” down the tree fragment, thus the name
“diving.”
Diving heuristics generally select one variable that should be integer-valued, for which the current solution is fractional. The heuristics then introduce a bound that forces the variable to be integer-valued, and solve the associated relaxed LP again. The method of choosing the variable to bound is the main difference between the diving heuristics. See Berthold [5], Section 3.1.
Randomized Rounding, RINS, and RENS
To find new integer-feasible points, intlinprog searches
the neighborhood of the current, best integer-feasible solution point (if
available) to find a new and better solution. See Danna, Rothberg, and Le Pape
[8]. Similarly, to find new integer-feasible points,
intlinprog takes the LP solution to the relaxed problem
at a node, and rounds the integer components in a way that attempts to maintain
feasibility. By taking randomized rounding steps,
intlinprog can sometimes find a new feasible point.
RENS, which stands for Relaxation Enforced Neighborhood Search, is another
search technique for finding integer-feasible points. See Berthold [5].
Branch-and-Bound
The branch-and-bound method constructs a sequence of subproblems that attempt to converge to a solution of the MILP. The subproblems give a sequence of upper and lower bounds on the solution fTx. These bounds are called the primal and dual bounds. For a minimization problem, the first upper bound (primal) is any feasible solution, and the first lower bound (dual) is the solution to the relaxed problem. For a maximization problem, the primal bound is the lower bound and the dual bound is the upper bound. Any solution to the linear programming relaxed problem has a lower objective function value than the solution to the MILP. Also, any feasible point xfeas satisfies
| fTxfeas ≥ fTx, | (1) |
because fTx is the minimum among all feasible points.
In this context, a node is an LP with the same objective function, bounds, and linear constraints as the original problem, but without integer constraints, and with particular changes to the linear constraints or bounds. The root node is the original problem with no integer constraints and no changes to the linear constraints or bounds, meaning the root node is the initial relaxed LP.
From the starting bounds, the branch-and-bound method constructs new
subproblems by branching from the root node. The branching step is taken
heuristically, according to one of several rules. Each rule is based on the idea
of splitting a problem by restricting one variable to be less than or equal to
an integer J, or greater than or equal to J+1. These two subproblems arise when
an entry in xLP, corresponding to an
integer specified in intcon, is not an integer. Here,
xLP is the solution to a relaxed
problem. Take J as the floor of the variable (rounded down), and J+1 as the
ceiling (rounded up). The resulting two problems have solutions that are larger
than or equal to
fTxLP,
because they have more restrictions. Therefore, for a minimization problem this
procedure potentially raises the lower bound.
After the algorithm branches, there are two new nodes to explore.
intlinprog skips the analysis of some subproblems by
comparing their objective function values with the current global bounds.
The branch-and-bound procedure continues, systematically generating subproblems to analyze and discarding the ones that will not improve an upper or lower bound on the objective, until one of these stopping criteria is met:
The problem is proved to be infeasible.
The objective function value reaches the
ObjectiveCutOfflimit.The algorithm exceeds the
MaxTimeoption.The difference between the lower and upper bounds on the objective function is less than the
AbsoluteGapToleranceorRelativeGapTolerancetolerances.The number of explored nodes exceeds the
MaxNodesoption.
For background about the branch-and-bound procedure, see Nemhauser and Wolsey [11] and Wolsey [13].
Iterative Display
When you select iterative display by setting the Display
option to the default "iter", the solver displays some of its
steps. HiGHS iterative display is more extensive and complicated than the
iterative display of other solvers. Furthermore, the HiGHS algorithm can restart
its branch-and-bound search, leading to an iterative display that also
restarts.
To select iterative display:
options = optimoptions("intlinprog",Display="iter"); [x,fval,exitflag,output] = intlinprog(f,intcon,A,b,Aeq,beq,... lb,ub,options)
Preamble. The iterative display begins by displaying the results of "presolve." The presolve algorithm reduces the complexity of the original problem by identifying and removing redundant rows and columns in the linear constraint matrices and performing related simplifications of the problem. For example,
Presolving model 18018 rows, 26027 cols, 248579 nonzeros 15092 rows, 24343 cols, 217277 nonzeros Objective function is integral with scale 1
Root Node Evaluation. The root node is the linear programming solution of the problem, not
considering any integer constraints. For a minimization problem, the
objective function value of the root node is a lower bound on the objective
function value of the solution to the problem including integer constraints.
For a minimization problem, the upper bound (if any) comes from a feasible
point with respect to all constraints. If there is no feasible point yet
found, the upper bound is Inf.
The iterative display shows the size of the problem after presolve.
Solving MIP model with: 15092 rows 24343 cols (24343 binary, 0 integer, 0 implied int., 0 continuous) 217277 nonzeros
binaryis the number of binary variables.integeris the number of integer variables.implied intis the number of variables that are implied to be integer. For example, ifx(1)is integer andx(1) + x(2) = 5, thenx(2)is implied to be integer.continuousis the number of continuous variables.
The total number of variables is the number of columns in the model.
Dynamic Constraint Creation. The software starts by creating "dynamic constraints," which have three headers in the iterative display:
Cuts — Number of active cuts
inLp — Number of non-model rows in the LP matrix
Confl. — Number of conflicts
The software can further extend the constraints by restarting the dynamic constraint creation, which can create more constraints by starting the creation process from a new state.
In the last step before beginning the branch-and-bound process, the software reports the results of symmetry detection and the number of generators and orbitopes found. For example, this is the portion of the iterative display that appears in the preamble and dynamic constraint creation phase:
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 1.5s
0 0 0 0.00% 0 inf inf 0 0 5 4934 3.6s
…
0 0 0 0.00% 0 inf inf 4630 553 289 140296 159.6s
0.0% inactive integer columns, restarting
Model after restart has 15090 rows, 24338 cols (24338 bin., 0 int., 0 impl., 0 cont.), and 217075 nonzeros
0 0 0 0.00% 0 inf inf 550 0 0 140296 160.5s
…
0 0 0 0.00% 0 inf inf 5602 524 260 277323 318.0s
6.2% inactive integer columns, restarting
Model after restart has 14185 rows, 22816 cols (22816 bin., 0 int., 0 impl., 0 cont.), and 200423 nonzeros
0 0 0 0.00% 0 inf inf 524 0 0 277323 318.6s
…
0 0 0 0.00% 0 inf inf 4683 535 525 408393 505.8s
Symmetry detection completed in 3.1s
Found 215 generators and 12 full orbitope(s) acting on 730 columnsFor information about symmetry detection and orbitopes, see Hojny, Pfetsch, and Schmitt [7] and Pfetsch and Rehn [12]. The software continues to add dynamic constraints as it proceeds with branch-and-bound steps, described next.
Branch-and-Bound. The branch-and-bound method constructs a sequence of subproblems that attempt to converge to a solution of the MILP. The subproblems give a sequence of upper and lower bounds on the solution fTx. For a minimization problem, the first upper bound is any feasible solution, and the first lower bound is the solution to the relaxed problem.
During the branch-and-bound procedure, intlinprog
gives the following iterative display.
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
72 0 2 1.56% 0 inf inf 4695 535 738 667009 686.9s
…
T 271 107 51 1.56% 0 449 100.00% 6051 271 1767 792786 776.8s
T 279 107 53 1.56% 0 439 100.00% 6053 271 1794 793241 777.5s
…
L 1223 538 295 1.98% 0 434 100.00% 6984 243 7628 1689k 1580.6s
…
1321 539 333 1.98% 0 434 100.00% 7029 243 8628 1898k 1650.6s
Restarting search from the root node
Model after restart has 13947 rows, 22426 cols (22426 bin., 0 int., 0 impl., 0 cont.), and 194633 nonzeros
1323 0 0 0.00% 0 434 100.00% 243 0 0 1902k 1653.5s
1323 0 0 0.00% 0 434 100.00% 243 76 10 1905k 1655.7s
…
1694 173 52 0.00% 0 434 100.00% 9411 318 2220 3205k 2584.3s
B 1710 167 55 0.00% 0 433 100.00% 9415 318 2263 3207k 2586.2s
1726 224 56 0.00% 0 433 100.00% 9999 378 2307 3237k 2608.7sThe leftmost column shows a code indicating how the new feasible point was found for that row of the display:
L— While solving a sub-MIP problem during primal heuristicsT— During tree search, while evaluating a nodeB— During branchingH— By heuristicsP— During startup, before solving MIPC— By central roundingR— By randomized roundingS— While solving an LPF— By Feasibility PumpU— From an unbounded LP
Finish. The reason that the algorithm stopped and a summary of results appear at the end of the iterative display.
Solving report
Status Time limit reached
Primal bound 14
Dual bound 0
Gap 100% (tolerance: 0.01%)
Solution status feasible
14 (objective)
0 (bound viol.)
1.7763568394e-15 (int. viol.)
0 (row viol.)
Timing 7200.02 (total)
3.08 (presolve)
0.00 (postsolve)
Nodes 5830
LP iterations 9963775 (total)
2657448 (strong br.)
324908 (separation)
2140011 (heuristics)
Solver stopped prematurely. Integer feasible point found.
Intlinprog stopped because it exceeded the time limit, options.MaxTime = 7200 . The intcon variables are integer within tolerance,
options.ConstraintTolerance = 1e-06.Status— Reason the iterations stoppedPrimal bound— Upper bound on objective function value for a minimization problem, lower bound for a maximization problemDual bound— Lower bound on objective function value for a maximization problem, upper bound for a minimization problemGap— Relative gap between primal and dual boundsSolution statusStatus of the solution
Objective function value, labeled
(objective)Maximum violation of the solution with respect to variable bounds, labeled
(bound viol.)Maximum violation of the integer-type variables from integer values, labeled
(int. viol.)Maximum violation of the solution with respect to the linear constraints, labeled
(row viol.)
Timing— Timing of various solution phases in secondsNodes— Number of nodes exploredLP iterations— Number of linear program iterations by phase and in total
References
[1] Achterberg, T.,Bixby, R.,
Gu, Z., Rothberg, E., and Weninger, D. Presolve Reductions in Mixed
Integer Programming. INFORMS J. Computing 32(2), 2019, pp. 473–506.
Available at https://doi.org/10.1287/ijoc.2018.0857.
[2] Achterberg, T., T. Koch
and A. Martin. Branching rules revisited. Operations Research
Letters 33, 2005, pp. 42–54. Available at https://opus4.kobv.de/opus4-zib/files/788/ZR-04-13.pdf.
[3] Andersen, E. D., and Andersen, K. D. Presolving in linear programming. Mathematical Programming 71, pp. 221–245, 1995.
[4] Atamtürk, A., G. L. Nemhauser, M. W. P. Savelsbergh. Conflict graphs in solving integer programming problems. European Journal of Operational Research 121, 2000, pp. 40–55.
[5] Berthold, T. Primal Heuristics for Mixed
Integer Programs. Technischen Universität Berlin, September 2006.
Available at https://opus4.kobv.de/opus4-zib/files/1029/Berthold_Primal_Heuristics_For_Mixed_Integer_Programs.pdf.
[6] Cornuéjols, G. Valid inequalities for mixed integer linear programs. Mathematical Programming B, Vol. 112, pp. 3–44, 2008.
[7] Hojny, C., Pfetsch, M. E.,
Schmitt, A. Extended formulations for column-constrained
orbitopes. TU Darmstadt, Department of Mathematics, Research Group
Optimization, Dolivostr. 15, 64293 Darmstadt, June 2017. Available at https://optimization-online.org/2017/06/6092/.
[8] Danna, E., Rothberg, E., Le Pape, C. Exploring relaxation induced neighborhoods to improve MIP solutions. Mathematical Programming, Vol. 102, issue 1, pp. 71–90, 2005.
[9] Hendel, G. New Rounding and Propagation Heuristics for Mixed Integer Programming. Bachelor's thesis at Technische Universität Berlin, 2011. PDF available at https://opus4.kobv.de/opus4-zib/files/1332/bachelor_thesis_main.pdf.
[10] Mészáros C., and Suhl, U. H. Advanced preprocessing techniques for linear and quadratic programming. OR Spectrum, 25(4), pp. 575–595, 2003.
[11] Nemhauser, G. L. and Wolsey, L. A. Integer and Combinatorial Optimization. Wiley-Interscience, New York, 1999.
[12] Pfetsch, M. E. and Rehn,
T. A computational comparison of symmetry handling methods for mixed
integer programs. Mathematical Programming Computation (2019)
11:37–93. Available at https://optimization-online.org/wp-content/uploads/2015/11/5209.pdf.
[13] Wolsey, L. A. Integer Programming. Wiley-Interscience, New York, 1998.
