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.
To enter a problem symbolically:
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.
constr:=[x<=1,y>=2,x+3y-z=2,3x-y+z<=8,-x+y<=5]; |
lpsolve(2x+y-z+4,constr) |
⎡ ⎣ | −4, | ⎡ ⎣ | x=0,y=5,z=13 | ⎤ ⎦ | ⎤ ⎦ |
⎡ ⎣ | −2, | ⎡ ⎣ | x=0,y=5,z=−4 | ⎤ ⎦ | ⎤ ⎦ |
⎡ ⎢ ⎢ ⎣ | − |
| , | ⎡ ⎢ ⎢ ⎣ | x= |
| ,y= |
| ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
constr:=[5x-10y<=20,2z-3y=6,-x+3y<=3]; |
lpsolve(-6x+4y+z,constr,x=1..20,y=0..inf) |
⎡ ⎢ ⎢ ⎣ | − |
| , | ⎡ ⎢ ⎢ ⎣ | x=18,y=7,z= |
| ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
To enter a problem in matrix form:
If the problem does not contain equality constraints, Aeq and beq may be omitted. For a problem that does not contain inequality constraints, empty lists must be entered in place of A and b.
Examples.
c:=[-2,1];A:=[[-1,1],[1,1],[-1,0],[0,-1]];b:=[3,5,0,0]; |
lpsolve(c,[A,b]) |
⎡ ⎣ | −10, | ⎡ ⎣ | 5,0 | ⎤ ⎦ | ⎤ ⎦ |
c:=[-2,5,-3];bl:=[2,3,1];bu:=[6,10,7/2]; |
lpsolve(c,[],[bl,bu]) |
⎡ ⎢ ⎢ ⎣ | − |
| , | ⎡ ⎢ ⎢ ⎣ | 6,3, |
| ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
c:=[4,5];Aeq:=[[-1,3/2],[-3,2]];beq:=[2,3]; |
lpsolve(c,[[],[],Aeq,beq]) |
⎡ ⎢ ⎢ ⎣ |
| , | ⎡ ⎢ ⎢ ⎣ | − |
| , |
| ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
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.
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.
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.
To solve pure integer programming problems, in which all variables are integers, use the option assume=integer or assume=lp_integer.
Example.
Input:
Output:
⎡ ⎣ | −41, | ⎡ ⎣ | x=4,y=3 | ⎤ ⎦ | ⎤ ⎦ |
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.
Example.
Input:
Output:
⎡ ⎣ | 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.
Input:
Output:
Use the assume=lp_nonnegint or assume=nonnegint option to get nonnegative integer values.
Example.
Input:
Output:
⎡ ⎣ | 12, | ⎡ ⎣ | x=1,y=2 | ⎤ ⎦ | ⎤ ⎦ |
When specifying MIP problems in matrix form, the lists corresponding to the options lp_integervariables and lp_binaryvariables should contain indices of the variables, as in the following example.
Example.
Input:
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]) |
Output:
⎡ ⎢ ⎢ ⎣ |
| , | ⎡ ⎢ ⎢ ⎣ | 1, |
| ,−1 | ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
You can also specify a range of indices instead of a list when there is too much variables. Example: lp_binaryvariables=0..99 means that all variables xi such that 0≤ i≤ 99 are binary.
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. Since the bounds of the problem variables are the only parameters that change accross the subproblems, the number of constraints does not change, which is a benefit of the upper-bounding technique built in the simplex algorithm. This makes the algorithm space-efficient.
A node-selection strategy can be chosen by using the lp_nodeselect option. Possible values are:
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.
A variable-selection strategy can be chosen by using the lp_varselect option. Possible values are:
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= |
| , b= |
| , c= |
| . |
The optimal solution is z∗=1854, x∗=(24,15,19,11,3,99,4,226). In the following table, different strategies are compared according to the number of examined nodes and total solving time (in seconds).
Node selection | Variable selection | Nodes examined | Time |
Best local bound | Last fractional | 13102 | 5.9 |
Best projection | Pseudocost | 23528 | 14.2 |
Hybrid | Most fractional | 43410 | 21.7 |
Breadth-first | Pseudocost | 77034 | 100.8 |
Note that the above comparison is problem-specific; it does not mean that lp_bestlocalbound+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.
Strong Gomory mixed integer cuts are generated at every node of the branch-and-bound tree and used to improve the objective value bound. After solving the relaxed subproblem by using the simplex method, at most one strong cut is generated and added to the subproblem which is subsequently reoptimized. This procedure is repeated until no more cuts can be generated.
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 cuts that is allowed be appended to a subproblem, you can 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.
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.
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 Gomory mixed integer cuts.
In the following example, we solve the MIP given above. The solver shows all progress and summary messages.
Example.
Input:
Output:
Constraint matrix has 4 rows, 8 columns, and 36 nonzeros |
Preprocessing... |
Constraint matrix has 4 rows, 8 columns, and 36 nonzeros |
Optimizing... |
Applying branch & bound method... |
76: Switched to pseudocost based branching |
8906: 909 nodes active, bound: 1574.31 |
9999: Incumbent solution found |
17052: 598 nodes active, bound: 1706.09, gap: 7.9779% |
23528: Tree is empty |
Summary: |
* 23528 subproblem(s) examined |
* max. tree size: 1021 nodes |
* 398 Gomory cut(s) applied |
⎡ ⎣ | 1854, | ⎡ ⎣ | 24,15,19,11,3,99,4,226 | ⎤ ⎦ | ⎤ ⎦ |
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 sophisticated LP/MIP solvers in floating-point arithmetic, designed to be very fast and to handle large 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_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. 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.
For example, we solve the following LP problem using the default settings.
subject to
|
Input:
lpsolve(1.06x1+0.56x2+3x3, |
[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05]) |
Output:
⎡ ⎣ | 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.
Input:
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]) |
Output:
⎡ ⎣ | 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.
Input:
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) |
Output:
⎡ ⎣ | 2255937.50731, | ⎡ ⎣ | x1=688490.256652,x2=2716245.85608,x3=1680.05195065 | ⎤ ⎦ | ⎤ ⎦ |
Linear (integer) programming problems can be loaded from MPS1 or CPLEX LP2 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:
Input:
Output:
⎡ ⎣ | 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.
Input:
Output:
⎡ ⎢ ⎢ ⎣ |
| , | ⎡ ⎢ ⎢ ⎣ | x1=0,x2= |
| ,x3=16 | ⎤ ⎥ ⎥ ⎦ | ⎤ ⎥ ⎥ ⎦ |
Next we force all variables to integral values.
Input:
Output:
⎡ ⎣ | 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 many 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 6.52.2).