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:
- 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
| x | ≤ 1, | | | | | | | | | |
y | ≥ 2, | | | | | | | | | |
x+3y−z | =2, | | | | | | | | | |
2x−y+z | ≤ 8, | | | | | | | | | |
−x+y | ≤ 5.
| | | | | | | | | |
|
lpsolve(2x+y-z+4,[x<=1,y>=2,x+3y-z=2,3x-y+z<=8,-x+y<=5]) |
|
| ⎡
⎣ | −4, | ⎡
⎣ | x=0,y=5,z=13 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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.
lpsolve(x+2y+3z,[x+y=1..5,y+z+1=2..4,x>=0,y>=0]) |
|
| ⎡
⎣ | −2, | ⎡
⎣ | x=0,y=5,z=−4 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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.
lpsolve(-x-y,[y<=3x+1/2,y<=-5x+2],assume=lp_nonnegative) |
|
| ⎡
⎢
⎢
⎣ | − | | , | ⎡
⎢
⎢
⎣ | x= | | ,y= | | ⎤
⎥
⎥
⎦ | ⎤
⎥
⎥
⎦ |
| | | | | | | | | | |
|
Bounds can be added separately for some variables. They should be
entered after the constraints. For example, minimize −6x+4y+z subject to
| 5x−10y | ≤ 20, | | | | | | | | | |
2z−3y | =6, | | | | | | | | | |
−x+3y | ≤ 3,
| | | | | | | | | |
|
where 1≤ x≤ 20 and y≥ 0.
lpsolve(-6x+4y+z,[5x-10y<=20,2z-3y=6,-x+3y<=3],x=1..20,y=0..inf) |
|
| ⎡
⎢
⎢
⎣ | − | | , | ⎡
⎢
⎢
⎣ | x=18,y=7,z= | | ⎤
⎥
⎥
⎦ | ⎤
⎥
⎥
⎦ |
| | | | | | | | | | |
|
Solving an LP problem in matrix form
To enter a problem in matrix form:
-
lpsolve takes
- 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.
c:=[-2,1]:; A:=[[-1,1],[1,1],[-1,0],[0,-1]]:; b:=[3,5,0,0]:;
lpsolve(c,[A,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.
lpsolve([-2,5,-3],[],[[2,3,1],[6,10,7/2]]) |
|
| ⎡
⎢
⎢
⎣ | − | | , | ⎡
⎢
⎢
⎣ | 6,3, | | ⎤
⎥
⎥
⎦ | ⎤
⎥
⎥
⎦ |
| | | | | | | | | | |
|
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.
lpsolve([4,5],[[],[],[[-1,3/2],[-3,2]],[2,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:
lpsolve(-5x-7y,[7x+y<=35,-x+3y<=6],assume=integer) |
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:
lpsolve(8x1+11x2+6x3+4x4,[5x1+7x2+4x3+3x4<=14],assume=lp_binary,maximize) |
|
| ⎡
⎣ | 21, | ⎡
⎣ | x1=0,x2=1,x3=1,x4=1 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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:
lpsolve(x+3y+3z,[x+3y+2z<=7,2x+2y+z<=11],
assume=lp_nonnegative,lp_maximize,lp_integervariables=[x,z] |
Use the assume=lp_nonnegint or assume=nonnegint
option to get nonnegative integer values. For example:
lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint) |
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:
c:=[2,-3,-5]:; A:=[[-5,4,-5],[2,5,7],[2,-3,4]]:; b:=[3,1,-2]:;
lpsolve(c,[A,b],lp_integervariables=[0,2]) |
|
| ⎡
⎢
⎢
⎣ | | , | ⎡
⎢
⎢
⎣ | 1, | | ,−1 | ⎤
⎥
⎥
⎦ | ⎤
⎥
⎥
⎦ |
| | | | | | | | | | |
|
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 xi 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∈ℤ+8 and
A= | ⎡
⎢
⎢
⎢
⎣ | 22 | 13 | 26 | 33 | 21 | 3 | 14 | 26 |
39 | 16 | 22 | 28 | 26 | 30 | 23 | 24 |
18 | 14 | 29 | 27 | 30 | 38 | 26 | 26 |
41 | 26 | 28 | 36 | 18 | 38 | 16 | 26
|
| ⎤
⎥
⎥
⎥
⎦ |
| ,
b= | | ,
c=[2,10,13,17,7,5,7,3]T.
|
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
n obj: z
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
n: m nodes active, bound: b⟨ gap: g⟩
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.
A:=[[22,13,26,33,21,3,14,26],[39,16,22,28,26,30,23,24],
[18,14,29,27,30,38,26,26],[41,26,28,36,18,38,16,26]]:;
b:=[7872,10466,11322,12058]:; c:=[2,10,13,17,7,5,7,3]:;
lpsolve(c,[[],[],A,b],assume=nonnegint,lp_verbose) |
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
|
| ⎡
⎣ | 1854, | ⎡
⎣ | 24,15,19,11,3,99,4,226 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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.
Minimize 1.06 x1+0.56 x2+3.0 x3
|
subject to
| 1.06 x1+0.015 x3 | ≥ 729824.87 | | | | | | | | | |
0.56 x2+0.649 x3 | ≥ 1522188.03 | | | | | | | | | |
x3 | ≥ 1680.05
| | | | | | | | | |
|
lpsolve(1.06x1+0.56x2+3x3,
[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05]) |
|
| ⎡
⎣ | 2255937.4968, | ⎡
⎣ | x1=688490.254009,x2=2716245.85277,x3=1680.05 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
If the requirement xk∈ℤ for k=1,2, the following result is
obtained from the MIP solver in GLPK.
lpsolve(1.06x1+0.56x2+3x3,
[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05],
lp_integervariables=[0,1]) |
|
| ⎡
⎣ | 2255938.37, | ⎡
⎣ | x1=688491,x2=2716246,x3=1680.05 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
The solution of the original problem can also be obtained with
the interior-point solver by using the
lp_method=lp_interiorpoint option.
lpsolve(1.06x1+0.56x2+3x3,
[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05],
lp_method=lp_interiorpoint) |
|
| ⎡
⎣ | 2255937.50731, | ⎡
⎣ | x1=688490.256652,x2=2716245.85608,x3=1680.05195065 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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:
Maximize
obj: 2x1+4x2+3x3
Subject To
c1: 3x1+5x2+2x3 <= 60
c2: 2x1+x2+2x3 <= 40
c3: x1+3x2+2x3 <= 80
General
x1 x3
End |
To find an optimal solution to the linear program specified in this file, enter:
|
| ⎡
⎣ | 71.8, | ⎡
⎣ | x1=0,x2=5.2,x3=17 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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.
lpsolve("somefile.lp",x2=5.4..inf,lp_method=exact) |
|
| ⎡
⎢
⎢
⎣ | | , | ⎡
⎢
⎢
⎣ | x1=0,x2= | | ,x3=16 | ⎤
⎥
⎥
⎦ | ⎤
⎥
⎥
⎦ |
| | | | | | | | | | |
|
Next we force all variables to integral values.
lpsolve("somefile.lp",x2=5.4..inf,assume=integer) |
|
| ⎡
⎣ | 69.0, | ⎡
⎣ | x1=0,x2=6,x3=15 | ⎤
⎦ | ⎤
⎦ |
| | | | | | | | | | |
|
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).