16.1.2  Solving (mixed integer) linear programming problems

The lpsolve command can solve linear programming problems, where a multivariate linear function needs to be maximized or minimized subject to linear (in)equality constraints, as well as (mixed) integer programming problems. You can either enter a problem directly (in symbolic or matrix form) or load it from a file in LP or (compressed) MPS format.

Solving an LP problem in symbolic form

lpsolve can solve LP problems in symbolic form in which the variables appear as symbols.

• lpsolve takes one mandatory argument and three optional arguments:
  • obj, a symbolic expression representing the objective function or a path to a file containing the LP problem.
  • Optionally, constr, list of linear constraints which may be equalities or inequalities or doubly bounded expressions in form a≤ f(x)≤ b entered as f(x)=a..b. (If obj is a file name, this option is omitted.)
  • Optionally, bd, a sequence of expressions of type x=a..b, specifying that the variable x is bounded below by a∈ℝ∪{−∞} and above by b∈ℝ∪{+∞}. If the bounds are stored in a two-column matrix, consider using the conversion routine box_constraints (see Section 6.6.9).
  • Optionally, opts, a sequence of solver settings in form option=value, where option may be one of:
    • lp_assume (or assume), which specifies a global constraint on variables and whose value can be one of:
      • lp_nonnegative (all variables are nonnegative)
      • lp_integer or integer (all variables are integers)
      • lp_binary (all variables are binary, i.e. 0 or 1)
      • lp_nonnegint or nonnegint (all variables are nonnegative integers)
      (by default, unset.)
    • lp_integervariables, whose value should be a list of identifiers or indices of integer variables (empty by default).
    • lp_binaryvariables, whose value should be a list of identifiers or indices of binary variables (empty by default).
    • lp_maximize (or maximize), whose value can be true or false setting the objective direction (by default false, meaning that the objective is minimized). You can enter only maximize, which is equivalent to maximize=true.
    • lp_method, setting the solver type, which can be one of: exact, float, lp_simplex (the default) or lp_interiorpoint.
    • lp_depthlimit, whose value can be a positive integer, which sets the maximum depth of the branch-and-bound tree (by default, unlimited).
    • lp_nodelimit, whose value can be a positive integer, which sets the maximum number of nodes in the branch-and-bound tree (by default, unlimited).
    • lp_iterationlimit, whose value can be a positive integer setting the maximum iterations of the simplex algorithm (by default, unlimited).

      If the maximum number of iterations is reached, the current feasible solution (not necessarily an optimal one) is returned.

    • lp_timelimit, whose value can be a positive real number, setting the maximum solving time in milliseconds (by default, unlimited).
    • acyclic, whose value can be either true or false, enabling/disabling the Bland rule safeguarding (by default, true).
    • lp_maxcuts, whose value can be a nonnegative integer setting the maximum GMI cuts per node (by default 5).
    • lp_gaptolerance, whose value can be a positive number, setting the relative integrality gap threshold (by default, 0).
    • lp_presolve, whose value can be either true or false, enabling/disabling the preprocessing step (by default, true).
    • lp_heuristic, whose value can be either true or false, enabling/disabling the rounding heuristic (by default, true).
    • lp_nodeselect, which sets the branching node selection strategy and whose value can be one of: lp_depthfirst, lp_breadthfirst, lp_bestlocalbound (the default), lp_hybrid or lp_bestprojection.
    • lp_varselect, which sets the branching variable selection strategy, whose value can be one of: lp_firstfractional, lp_lastfractional, lp_mostfractional or lp_pseudocost (the default).
    • lp_verbose, whose value can be true or false (by default false). You can enter only lp_verbose, which is equivalent to lp_verbose=true.
• lpsolve(obj ⟨,constr,bd,opts ⟩) returns a list [optimum,soln], where optimum is the minimum/maximum value of the objective function and soln is the list of coordinates corresponding to the point at which the optimal value is attained, i.e. the optimal solution. If there is no feasible solution, an empty list is returned. When the objective function is unbounded, optimum is returned as +infinity (for maximization problems) or -infinity (for minimization problems). If an error is experienced while solving, then the process is terminated and undef is returned.

The given objective function is minimized by default. To maximize it, include the option lp_maximize=true or lp_maximize or simply maximize. Also note that all variables are, unless specified otherwise, assumed to be continuous and unrestricted in sign.

Examples

Minimize 2x+y−z+4 subject to

Therefore, the minimum value of f(x,y,z)=2 x+y−z+4 is equal to −4 under the given constraints. The optimal value is attained at point (x,y,z)=(0,5,13).

Constraints may also take the form expr=a..b for bounded linear expressions. For example, minimize x+2y+3z subject to 1≤ x+y≤ 5 and 2≤ y+z+1≤ 4, where x,y≥ 0.

Use the assume=lp_nonnegative option to specify that all variables are nonnegative. It is easier than entering the nonnegativity constraints explicitly. For example, minimize −x−y subject to y≤ 3x+1/2 and y≤ −5x+2, where x,y≥ 0.

Bounds can be added separately for some variables. They should be entered after the constraints. For example, minimize −6x+4y+z subject to

where 1≤ x≤ 20 and y≥ 0.

Solving an LP problem in matrix form

To enter a problem in matrix form:

• lpsolve takes
  • obj, a vector of coefficients representing the objective function.
  • constr, a list [A,b,A_eq,b_eq] such that objective function obj^T· x is to be minimized/maximized subject to constraints Ax ≤ b and A_eqx=b_eq.

    If the problem does not contain equality constraints, A_eq and b_eq may be omitted. For a problem that does not contain inequality constraints, empty lists must be entered in place of A and b.

  • bd, a list of two vectors [b_l,b_u] of the same length as c such that b_l≤ x ≤ b_u. These vectors may contain +infinity or -infinity.
  • opts, as before.
• lpsolve(obj ⟨,constr,bd,opts ⟩) returns a list [optimum,soln] as before. optimum is the minimum/maximum value of the objective function and soln is the list of coordinates corresponding to the point at which the optimal value is attained, i.e. the optimal solution. If there is no feasible solution, an empty list is returned. When the objective function is unbounded, optimum is returned as +infinity (for maximization problems) or -infinity (for minimization problems). If an error is experienced while solving (terminating the process), undef is returned.

Examples

We minimize −2x+y subject to −x+y≤ 3, x+y≤ 5, x≥ 0 and y≥ 0 by providing the parameters A and b.

You can enter variable bounds as the third argument. The following example minimizes −2x+5y−3z subject to 2≤ x≤ 6, 3≤ y≤ 10 and 1≤ z≤ 7/2.

To solve problems with only equality constraints, enter empty lists in places of A and b. The following example minimizes 4x+5y subject to −x+3/2 y=2 and −3x+2y=3.

Simplex method implementation.

Only the two-phase primal simplex method is implemented for lpsolve and it uses the upper-bounding technique. Basic columns are not kept in the simplex tableau, which therefore contains only the columns of non-basic variables. To prevent cycling, an adaptation of Bland’s rule is used. A preprocessing routine, helping to reduce the size of the problem, is available and enabled by default (you can disable it by typing lp_presolve=false). All computation is done by using arbitrary-precision integer arithmetic, which is always exact but slower than the floating-point arithmetic. Note that problem data should be rational. (To solve LP problems in floating-point arithmetic, see Section 16.1.2.)

Cycling in simplex algorithm may happen, although it is rare in practice. Bland rule safeguarding is used by default in order to prevent cycling. However, Bland’s pivoting rule is known to converge slowly; therefore it may slow down the computation in problems which would otherwise not cycle. To disable the safeguarding, use the option acyclic=false.

Solving MIP (Mixed Integer Programming) problems

The lpsolve command allows restricting (some) variables to integer values. Such problems, called (mixed) integer programming problems, are solved by applying the branch-and-bound method.

Examples

To solve pure integer programming problems, in which all variables are integers, use the option assume=integer or assume=lp_integer. For example:

Use the option assume=lp_binary to specify that all variables are binary, i.e. the only allowed values are 0 and 1. These usually represent false and true, respectively, giving the variable a certain meaning in a logical context. For example:

To solve mixed integer problems, where some variables are integers and some are continuous, use the option keywords lp_integervariables to specify integer variables and/or lp_binaryvariables to specify binary variables. For example:

Use the assume=lp_nonnegint or assume=nonnegint option to get nonnegative integer values. For example:

When specifying MIP problems in matrix form, the lists corresponding to the options lp_integervariables and lp_binaryvariables should contain indices of the variables. For example:

You can also specify a range of indices instead of a list when there is too much variables. For example, lp_binaryvariables=0..99 means that all variables x_i such that 0≤ i≤ 99 are binary.

Branch-and-bound implementation.

The branch-and-bound algorithm generates a binary tree of subproblems, called nodes, by branching on integer variables with fractional values. Leaf nodes of the tree, called active nodes, are stored in a list. In each iteration of the algorithm, an active node is selected, branched on a particular variable into two new nodes, and subsequently removed from the list. A node in which no branching is possible represents a feasible solution. The corresponding objective value is used to fathom nodes which cannot possibly lead to a better solution. The algorithm terminates when there is no space left for improvement.

If presolving is enabled, then basic preprocessing is done at each node of the tree, except when lp_presolve=root is set, in which case only the root node is processed. Additionally, after a non-integer-feasible solution with better objective value than the current incumbent is obtained by solving the linear relaxation, a rounding heuristic is applied in attempt to achieve integral feasibility. The heuristic is enabled by default; you can disable it by setting lp_heuristic to false.

Node-selection strategies.

A node-selection strategy can be chosen by using the lp_nodeselect option. Possible values are:

• lp_bestlocalbound, which chooses an active node having the best bound for the objective value,
• lp_depthfirst, which chooses the newest active node,
• lp_hybrid, which combines the above two strategies,
• lp_breadthfirst, which chooses the oldest active node,
• lp_bestprojection, which chooses a node with the best simple projection.

By default, the lp_bestlocalbound strategy is used. The lp_hybrid strategy works as follows: until an incumbent solution is found, the solver uses the lp_depthfirst strategy, “diving” into the tree as an incumbent solution is more likely to be located deeply. When an incumbent is found, the solver switches to the lp_bestlocalbound strategy in attempt to close the integrality gap as quickly as possible.

Variable-selection strategies.

A variable-selection strategy can be chosen by using the lp_varselect option. Possible values are:

• lp_firstfractional, which chooses the first fractional variable,
• lp_lastfractional, which chooses the last fractional variable,
• lp_mostfractional, which chooses a variable with fractional part closest to 0.5,
• lp_pseudocost, which chooses the variable which had the greatest impact on the objective value in previous branchings.

By default, the lp_pseudocost strategy is used. However, since a pseudocost-based choice cannot be made until all integer variables have been branched upon at least one time in each direction, the lp_mostfractional strategy is used until that condition is fulfilled.

Using an appropriate combination of node/variable selection strategies may significantly reduce the number of subproblems needed to be examined when solving a particular MIP problem, as illustrated in the following example.

Example

Minimize z=c·x subject to A x=b, where x∈ℤ₊⁸ and

The optimal solution is z^∗=1854, x^∗=(24,15,19,11,3,99,4,226).

In Table 16.1, different strategies are compared according to the number of examined nodes and total solving time (in seconds). Note that the above comparison is problem-specific; it does not mean that lp_bestlocalbound with lp_lastfractional is always the best strategy. Usually, one has to experiment with different combinations to find which one is optimal for the given problem. However, the strategies which use larger amounts of information generally perform better. Therefore, it is reasonable to assume that lp_bestlocalbound will be more appropriate than lp_breadthfirst, for instance.

---

Node selection	Variable selection	Nodes examined	Time
Best local bound	Last fractional	13102	4.8
Best projection	Most fractional	26238	11.7
Hybrid	First fractional	55046	19.5
Depth-first	Pseudocost	131466	36.2

Table 16.1: Comparison of different variable and subproblem selection strategies for lpsolve

---

Cutting planes.

Gomory mixed integer (GMI) cuts are generated at every node of the branch-and-bound tree to improve the objective value bound. After solving the relaxed subproblem using the simplex method, GMI cuts are added to the subproblem which is subsequently reoptimized. This procedure is repeated until no useful cuts can be generated or until lp_maxcuts limit is reached.

Simplex reoptimizations are fast because they start with the last feasible basis; however, applying cuts makes the simplex tableau larger, which may slow the computation down. To limit the number of GMI cuts that is allowed be appended to a subproblem, use lp_maxcuts option, setting it either to zero (which disables the cut generation altogether) or to some positive integer. You may set it to +infinity as well, thus allowing any number of cuts to be applied to a node. By default, lp_maxcuts equals to 5.

Stopping criteria.

There are several ways to force the branch-and-bound algorithm to stop prematurely when the execution takes too much time. You can set lp_timelimit to an integer which defines the maximum number of milliseconds allowed to find an optimal solution. Other ways are to set lp_nodelimit or lp_depthlimit to limit the number of nodes generated in the branch-and-bound tree or its depth, respectively. Finally, you can set lp_gaptolerance to some positive value, say t>0, which terminates the algorithm after finding an incumbent solution and proving that the corresponding objective value differs from the optimum value for less than t· 100 % . This is done by monitoring the size of the integrality gap, i.e. the difference between the current incumbent objective value and the best objective value bound among active nodes.

If the branch-and-bound algorithm terminates prematurely, a warning message indicating the cause is displayed. The incumbent solution, if any, is returned as the result, else the problem is declared to be infeasible.

Displaying detailed output.

Typing lp_verbose=true or lp_verbose when specifying options for lpsolve causes detailed messages to be printed during and after solving a LP problem. During the simplex algorithm, a status report in form

is printed every two seconds, where n is the number of simplex iterations and z is the current objective value. If the line is prefixed with the asterisk character (*), it means that the solver is optimizing the given objective (Phase 2); otherwise, the solver is searching for a feasible basis (Phase 1), in which case z is a relative value in percentages (when it reaches zero, Phase 2 is initiated). Note that values of z reported in Phase 2 may not correspond to the actual values if presolving is enabled.

If the problem contains integrality constraints, then the simplex algorithm messages are not printed. Instead, during the branch-and-bound algorithm, a status report in form

is displayed every 5 seconds, where n is the number of already examined subproblems, b is the lower resp. upper bound (for minimization resp. maximization), and g is the relative integrality gap. Also, a message is printed every time the incumbent solution is found or updated, as well as when the solver switches to pseudocost-based branching. After the algorithm is finished, a summary is displayed containing the total number of examined subproblems, the number of most nodes being active at the same time, and the number of applied GMI cuts along with the respective average objective value improvement.

Example

Here we solve the MIP given above. The solver shows all progress and summary messages.

Constraint matrix has 4 rows, 9 columns, and 36 nonzeros
Variables: 0 continuous, 8 integer (0 binary)
Constraints: 4 equalities, 0 inequalities
Optimizing...

Starting branch & bound...

11310: 147 nodes active, bound: 1793.43
12709: Incumbent solution found: 1854
20646: 931 nodes active, bound: 1841.81, gap: 0.657343%
23836: Tree is empty
Summary:
* 23836 subproblem(s) examined
* max. tree size: 1256 nodes
* 13 GMI cut(s) applied

Solving problems in floating-point arithmetic

The lpsolve command provides, in addition to its own exact solver implementing the primal simplex method with the upper-bounding technique, an interface to GLPK (GNU Linear Programming Kit) library which contains LP/MIP solvers operating in floating-point arithmetic, designed to be fast and to handle larger problems. Choosing between the available solvers is done by setting lp_method option.

By default, lp_method is set to lp_simplex, which solves the problem by using the native solver, but only if all problem coefficients are exact. If at least one of them is approximative (a floating-point number), then the GLPK solver is used instead.

Setting lp_method to exact forces the solver to perform exact computation even when some coefficients are inexact (they are converted to rational equivalents before applying the simplex method). In this case the native solver is used.

Setting lp_method to float forces lpsolve to use the GLPK simplex solver. If a (mixed) integer problem is given, then the branch-and-cut algorithm in GLPK is used. The parameters can be controlled by setting the lp_timelimit, lp_iterationlimit, lp_gaptolerance, lp_maxcuts, lp_heuristic, lp_verbose, lp_nodeselect, and lp_varselect options. If lp_varselect is not set, then the Driebeek–Tomlin heuristic is used, and if lp_nodeselect is not set, then the best-local-bound selection method is used. If lp_maxcuts is greater than zero, then GMI/MIR cut generation is enabled, else it is disabled. If the problem contains binary variables, then cover/clique cut generation is enabled, else it is disabled. If lp_heuristic=false, then the simple rounding heuristic in GLPK is disabled (by default it is enabled). Finally, lp_verbose=true enables GLPK messages, which are useful for monitoring solver’s progress.

Setting lp_method to lp_interiorpoint uses the primal-dual interior-point algorithm in GLPK. The only optional argument that affects this kind of solver is lp_verbose. The interior-point algorithm is faster than the simplex method for large and sparse LP problems. Note, however, that it does not support integrality constraints.

Example

We solve the following LP problem using the default settings.

subject to

If the requirement x_k∈ℤ for k=1,2, the following result is obtained from the MIP solver in GLPK.

The solution of the original problem can also be obtained with the interior-point solver by using the lp_method=lp_interiorpoint option.

Loading problems from files

Linear (integer) programming problems can be loaded from MPS or CPLEX LP format files. The file name should be a string passed as the obj parameter. If the file name has extension .lp, then CPLEX LP format is assumed. If the extension is .mps resp. .gz, then MPS resp. gzipped MPS format is assumed.

For example, assume that the file somefile.lp is stored in the working directory and that it contains the following lines of text:

To find an optimal solution to the linear program specified in this file, enter:

You can provide additional variable bounds, assumptions, and options alongside the file name, as in the examples below. Note that the original constraints (those which are read from file) cannot be removed.

Next we force all variables to integral values.

It is advisable to use only (capital) letters, digits and underscores when naming variables in an LP file, although the corresponding format allows more characters. That is because these names are converted to giac identifiers during the loading process.

Note that the coefficients of a problem loaded from file are always floating-point values. Therefore, to solve the problem with the native solver, you have to use the option lp_method=exact (see Section 16.1.2).