Solving general linear programming problems: lpsolve

6.52.2 Solving general linear programming problems: lpsolve

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 enter a problem directly (in symbolic or matrix form) load it from a file in LP or (gzipped) MPS format.

Solving an LP problem in symbolic form

To enter a problem symbolically:

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 bounded expressions entered as expr=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.
- Optionally, opts, a sequence of solver settings in form option=value, where option may be one of:
  - 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) (by default, lp_integervariables=[]).
  - lp_binaryvariables, whose value should be a list of identifiers or indices (of binary variables) (by default, lp_binaryvariables=[]).
  - lp_maximize (or maximize), whose value can be true or false setting the objective direction (by default, lp_maximize=false, meaning that the objective is minimized).
  - lp_method, setting the solver type, whose value can be one of:
    - exact
    - float
    - lp_simplex
    - lp_interiorpoint
    (by default, lp_method=lp_simplex).
  - 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 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 number, setting the maximum solving time in milliseconds (by default, unlimited).
  - lp_maxcuts, whose value can be a nonnegative integer setting the maximum GMI cuts per node, (by default, lp_maxcuts=5).
  - lp_gaptolerance, whose value can be a positive number, setting the relative integrality gap threshold (by default, lp_gaptolerance=0).
  - lp_nodeselect, which sets the branching node selection strategy and whose value can be one of:
    - lp_depthfirst
    - lp_breadthfirst
    - lp_hybrid
    - lp_bestprojection
    (by default, lp_nodeselect=lp_hybrid).
  - lp_varselect, which sets the branching variable selection strategy, whose value can be one of
    - lp_firstfractional
    - lp_lastfractional
    - lp_mostfractional
    - lp_pseudocost
    (by default, lp_varselect=lp_pseudocost).
  - lp_verbose, whose value can be true or false (by default: lp_verbose=false).
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 (terminating the process), 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.

Solve the problem specified in (11):
Input:

constr:=[x<=1,y>=2,x+3y-z=2,3x-y+z<=8,-x+y<=5];

lpsolve(2x+y-z+4,constr)

Output:
⎡
⎣ −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.
Input:
lpsolve(x+2y+3z,[x+y=1..5,y+z+1=2..4,x>=0,y>=0])
Output:
⎡
⎣ −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.
Input:
lpsolve(-x-y,[y<=3x+1/2,y<=-5x+2],assume=lp_nonnegative)
Output:
⎡
⎢
⎢
⎣ −
5

4

, ⎡
⎢
⎢
⎣ x=
3

16

,y=
17

16

⎤
⎥
⎥
⎦ ⎤
⎥
⎥
⎦
Bounds can be added separately for some variables. They should be entered after constraints.
Input:

constr:=[5x-10y<=20,2z-3y=6,-x+3y<=3];

lpsolve(-6x+4y+z,constr,x=1..20,y=0..inf)

Output:
⎡
⎢
⎢
⎣ −
133

2

, ⎡
⎢
⎢
⎣ x=18,y=7,z=
27

2

⎤
⎥
⎥
⎦ ⎤
⎥
⎥
⎦

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.

Input:

c:=[-2,1];A:=[[-1,1],[1,1],[-1,0],[0,-1]];b:=[3,5,0,0];

lpsolve(c,[A,b])

Output:
⎡
⎣ −10, ⎡
⎣ 5,0 ⎤
⎦ ⎤
⎦
Input:

c:=[-2,5,-3];bl:=[2,3,1];bu:=[6,10,7/2];

lpsolve(c,[],[bl,bu])

Output:
⎡
⎢
⎢
⎣ −
15

2

, ⎡
⎢
⎢
⎣ 6,3,
7

2

⎤
⎥
⎥
⎦ ⎤
⎥
⎥
⎦

Input:

c:=[4,5];Aeq:=[[-1,3/2],[-3,2]];beq:=[2,3];

lpsolve(c,[[],[],Aeq,beq])

Output:
⎡
⎢
⎢
⎣
26

5

, ⎡
⎢
⎢
⎣ −
1

5

,
6

5

⎤
⎥
⎥
⎦ ⎤
⎥
⎥
⎦

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.

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

Example.
Input:

lpsolve(-5x-7y,[7x+y<=35,-x+3y<=6],assume=integer)

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 logical context.

Example.
Input:

lpsolve(8x1+11x2+6x3+4x4,[5x1+7x2+4x3+3x4<=14],assume=lp_binary,maximize)

Output:

⎡
⎣

21,

⎡
⎣

x₁=0,x₂=1,x₃=1,x₄=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:

lpsolve(x+3y+3z,[x+3y+2z<=7,2x+2y+z<=11],
assume=lp_nonnegative,lp_maximize, lp_integervariables=[x,z])

Output:

[10,[x=1,y=0,z=3]]

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

Example.
Input:

lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint)

Output:

⎡
⎣

12,

⎡
⎣

x=1,y=2

⎤
⎦

When specifying MIP problems in matrix form, lists corresponding to the options lp_integervariables and lp_binaryvariables are populated with variable indices, like 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:

⎡
⎢
⎢
⎣

19,

⎡
⎢
⎢
⎣

,−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 x_i such that 0≤ i≤ 99 are binary.

About lpsolve

Implementation details.

The branch and bound algorithm by definition generates a binary tree of subproblems by branching on integer variables with fractional values. The lpsolve command features an implementation which stores only active nodes of a branch and bound tree in a list, thus saving a lot of space. Also, since variable bounds are the only parameters that change during the branch and bound algorithm, the number of constraints does not rise with depth, which is the benefit of the upper-bounding technique built in the simplex algorithm. Therefore a steady speed and minimal resource usage is always maintained, no matter how long the execution time is. This allows for solving problems that require tens or hundreds of thousands of nodes to be generated before finding an optimal solution.

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.

Branching strategies.

At every iteration of the branch and bound algorithm, a node must be selected for branching on some variable that has a fractional optimal value for the corresponding relaxed subproblem. There exist different methods for making such decisions, called branching strategies. Two types of branching strategies exist: the node selection and variable selection strategies.

The node selection strategy can be set by using the lp_nodeselect option. Possible values are:

lp_breadthfirst, which chooses the active node which provides the best bound for the objective value,
lp_depthfirst, which chooses the deepest active node and break ties by selecting the node providing the best bound,
lp_hybrid, which combines the above two strategies,
lp_bestprojection, which chooses the node with best simple projection.

By default, the lp_bestprojection strategy is used. Another sophisticated strategy is lp_hybrid: before 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_breadthfirst strategy to try to close the integrality gap as quickly as possible.

The Variable selection strategy can be set 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 the 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 pseudocost-based choice cannot be made before 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 the right combination of branching strategies may significantly reduce the number of subproblems needed to be examined when solving a particular MIP problem. However, what is “right” varies from problem to problem. Default strategies are the most sophisticated (as they use the available data most extensively) and usually the most effective ones. But that is not always the case, as illustrated by the following example.

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

⎡
⎢
⎢
⎢
⎣

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

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦

When using the default settings, about 24000 subproblems need to be examined before an optimal solution is found. When lp_nodeselect is set to lp_breadthfirst the solver needs to examine only about 20000 subproblems, but when set to lp_hybrid (a strategy which in general performs better) it examines about 111000 nodes in total.

Cutting planes.

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 with the simplex method, at most one strong cut is generated and added to the subproblem which is subsequently reoptimized. Simplex reoptimizations are fast because they start with the last feasible basis, but applying cuts makes the simplex tableau larger, hence applying many of them may actually slow the computation down. To limit the number of cuts that can be applied to a subproblem, you can use lp_maxcuts option, setting it either to zero (which disables cut generation altogether) or to some positive integer. Also, you may set it to +infinity, which means that any number of cuts may be applied to any node. By default, lp_maxcuts equals to 5.

Displaying detailed output.

By typing lp_verbose=true or simply lp_verbose when specifying options for lpsolve, detailed messages are printed during and after solving an MIP problem. During the branch and bound algorithm a status report in form

n: m nodes active, lower bound: lb ⟨ integrality gap: g⟩

is displayed every 5 seconds, where n is the number of already examined subproblems. Also, a report 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, i.e. when an optimal solution is found, 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, two nonnegative integers x₁ and x₂ are found such that 1867 x₁+1913 x₂=3618894 and x₁+x₂ is minimal. The solver shows all progress and summary messages.

Example.
Input:

lpsolve(x1+x2,[1867*x1+1913*x2=3618894],assume=nonnegint,lp_verbose=true)

Output:

Optimizing...

Applying branch&bound method to find integer feasible solutions...

3937: Incumbent solution found

Summary:

* 3938 subproblem(s) examined

* max. tree size: 1 nodes

* 0 Gomory cut(s) applied

⎡
⎣

1916,

⎡
⎣

x₁=1009,x₂=907

⎤
⎦

Solving problems in floating-point arithmetic

The lpsolve command provides, in addition to its own exact solver implementing the primal simplex method with upper-bounding technique, an interface to the 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 using primal simplex method, but performing exact computation only when all problem coefficients are exact. If at least one of them is approximative (a floating-point number), the GLPK solver is used instead (see below).

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).

Specifying lp_method=float forces lpsolve to use the floating-point solver. If a MIP problem is given, it is combined with the branch and cut algorithm. The GLPK simplex solver parameters can be controlled by setting the lp_timelimit, lp_gaptolerance and lp_varselect options. If the latter is not set, the Driebeek–Tomlin heuristic is used by default (see the GLPK manual for details). If lp_maxcuts is greater than zero, GMI and MIR cut generation is enabled, else it is disabled. If the problem contains binary variables, cover and clique cut generation is enabled, else it is disabled. Finally, lp_verbose=true enables detailed messages.

Setting lp_method to lp_interiorpoint uses the primal-dual interior-point algorithm which is part of GLPK. The only parameter that can be controlled via options is the verbosity level.

For example, try to solve the following LP problem using the default settings.

Minimize 1.06 x₁+0.56 x₂+3.0 x₃

subject to

1.06 x1+0.015 x3	≥ 729824.87
0.56 x₂+0.649 x₃	≥ 1522188.03
x₃	≥ 1680.05
x_k	≥ 0 for k=1,2,3

Input:

lpsolve(1.06x1+0.56x2+3x3,

[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05],

assume=lp_nonnegative)

Output:

⎡
⎣

2255937.4968,

⎡
⎣

x₁=688490.254009,x₂=2716245.85277,x₃=1680.05

⎤
⎦

If assume=nonnegint is used for the same problem, i.e. when x_k∈ℤ₊ for k=1,2,3, the following result is obtained by GLPK MIP solver: Input:

lpsolve(1.06x1+0.56x2+3x3,

[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05],

assume=nonnegint)

Output:

⎡
⎣

2255940.66,

⎡
⎣

x₁=688491.0,x₂=2716245.0,x₃=1681.0

⎤
⎦

The solution of the original problem can also be obtained with the interior-point solver by including lp_method=lp_interiorpoint after assume=lp_nonnegative: Input:

lpsolve(1.06x1+0.56x2+3x3,

[1.06x1+0.015x3>=729824.87,0.56x2+0.649x3>=1522188.03,x3>=1680.05],

assume=lp_nonnegative,lp_method=lp_interiorpoint)

Output:

⎡
⎣

2255937.50731,

⎡
⎣

x₁=688490.256652,x₂=2716245.85608,x₃=1680.05195065

⎤
⎦

Loading a problem from a file

Linear (integer) programming problems can be loaded from MPS or CPLEX LP format files (these formats are described in GLPK manual, Appendices B and C). The file name should be a string passed as the obj parameter. If the file name has extension “lp”, CPLEX LP format is assumed, and if the extension is “mps” or “gz”, MPS or gzipped MPS format is assumed.

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

\* Problem: short *\

Maximize
 obj: + 0.6 x1 + 0.5 x2

Subject To
 c1: + x1 + 2 x2 <= 1
 c2: + 3 x1 + x2 <= 2

Bounds
 x1 free
 x2 free

End

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

lpsolve("/path/to/file/somefile.lp")

Output:

⎡
⎢
⎢
⎣

x₁=

,x₂=

⎤
⎥
⎥
⎦

You can provide additional variable bounds and options alongside the file name. Note that the original constraints (those which are read from file) cannot be removed.
Input:

lpsolve("/path/to/file/somefile.lp",assume=integer)

Output:

⎡
⎣

x₁=1,x₂=−1

⎤
⎦

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.

Warning! Problems that are too large won’t be loaded. More precisely, those with n_v· n_c>10⁵, where n_v is the number of variables and n_c is the number of constraints, won’t be loaded.