Back

Software for solving linear programming problems

Linear programming problems and mixed-integer linear programming problems are rarely solved by hand even though you will be required to do so for some problems in this course so that you can become familiar with practical aspects of the theory. Problems from industrial applications often have thousands (and sometimes millions) of variables and constraints. Fortunately, there exist a number of commercial as well as open-source solvers that can handle such large-scale problem. We will now look at a number of options for solving LP problems using a computer.

Wolfram Alpha

You can solve linear programming problems using Wolfram Alpha. A bonus is that you can usually obtain exact rational solutions instead of floating-point ones. For example, to solve \[\begin{array}{rrcrcl} \min & 2x & + & y \\ \mbox{s.t.} & 3x & - & y & = & 2 \\ & x & + & 2y & \geq & 1 \\ & x & , & y & \geq & 0 \end{array}\] one can enter the following into Wolfram Alpha: 


Minimize[{2x+y, 3x - y == 2 && x + 2y >= 1 && x >= 0 && y >= 0}, {x,y}]

The result can be viewed here. It shows that \((x,y) = \left(\frac{5}{7},\frac{1}{7}\right)\) is an optimal solution.

Microsoft Excel or LibreOffice

If you have access to Microsoft Excel, you can also use it to solve all kinds of optimization problems. More information can be found here.

Since all students at Carleton have access to Excel, this will most likely be the first choice for solving larger LP problems using a computer in this course. Note that the solver add-in needs to be activated before it can be used.

An Excel spreadsheet that has the basic details set up for the problem \[\begin{array}{rrcrcrcl} \min & x & + & 2y & + & 4z\\ \mbox{s.t.} & 3x & - & y & - & z & = & 2 \\ & x & + & 2y & + & 2z & \geq & 1 \\ & x & , & y & , & z & \geq & 0 \end{array}\] can be found here.

If you do not have access to Microsoft Excel, you can use LibreOffice instead. There are slight differences between the two but either is sufficient for solving most of the problems that arise this course. Watch this video to see how to use LibreOffice to solve a linear programming problem.

Using JuMP in Julia

The recommended modelling tool for optimization problems in this course is the JuMP modeling language for the Julia Programming Language. You can download the latest version of Julia (v1.9.1 as of July 5, 2023) here. If you are a Mac user, consider installing via the package manager Homebrew which can save some headache (click here for instructions.)

There are two ways to run Julia. One is terminal based. The other is via Visual Studio Code or VSCodium with the Julia plugin. The latter is recommended unless you are a macOS/Linux user very familiar with terminal-based development or you are using it to launch the interactive environment Pluto. Instructions for installing Visual Code Studio for Julia programming can be found at Visual Studio Code and Julia in VS Code.

You can find a video on Julia development in VS Code here. Note that the video is missing the important fact that Ctrl-Shift-P (Command-Shift-P on macOS) brings up the command palette where you can type things like julia: start repl and julia: show plot.

Once everything is installed, install the JuMP and HiGHS packages via the following:

  1. Start Julia REPL (see this video or the VS Code video above if you are not sure)
  2. Run the following code:
    using Pkg
    Pkg.add("JuMP")
    Pkg.add("HiGHS")
You can solve the first example above using the following Julia code: 

using JuMP, HiGHS
model = Model();
set_optimizer(model, HiGHS.Optimizer)

@variable(model, x >= 0)
@variable(model, y >= 0)

@constraint(model, c1, 3x -  y == 2) 
@constraint(model, c2,  x + 2y >= 1)
@objective(model, Min, 2x + y)

optimize!(model) 
if termination_status(model) == MOI.OPTIMAL
  println("optval = ", objective_value(model))
  println("x = ", value(x))
  println("y = ", value(y))
else
  println("Unable to find an optimal solution.")
end

Using Google OR-Tools in Jupyter notebooks

There are also robust tools for solving mathematical optimization problems in Python. Here, we describe how to use Google OR-Tools.

First, you must install Jupyter. Instructions can be found here. Then install Google OR-Tools following the instructions.

Once installation is complete, you can launch Jupyter notebook. Say we want to solve the following: \[\begin{array}{rrcrcl} \max & x & - & 2y \\ \mbox{s.t.} & 3x & - & y & = & 2 \\ & x & + & 2y & \leq & 1 \\ & & & y & \in & \{0,1\} \end{array}\] We can enter the following in a cell and execute to solve the problem: 


from ortools.linear_solver import pywraplp

solver = pywraplp.Solver('simple_example',
                         pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

infinity = solver.infinity()

x = solver.NumVar(0., infinity, 'x')
y = solver.BoolVar('y')

solver.Add( 3.*x - y == 2. )
solver.Add( x + 2.*y <= 1.)

solver.Maximize( x - 2. * y )

result_status = solver.Solve()
if (result_status == pywraplp.Solver.OPTIMAL):
    print('Solution:')
    print('Objective value = ', solver.Objective().Value())
    print('x = ', x.solution_value())
    print('y = ', y.solution_value())
    
    # The following are other possible result statuses for your information
elif (result_status == pywraplp.Solver.INFEASIBLE):
    print("Problem is infeasible")
elif (result_status == pywraplp.Solver.UNBOUNDED):
    print("Problem is unbounded")
elif (result_status == pywraplp.Solver.ABNORMAL): 
    print("Problem is abnormal")
elif (result_status == pywraplp.Solver.NOT_SOLVED): 
    print("Problem is not solved")
Remarks

  1. The first argument to pywraplp.Solver can be replaced with a name of your choice specified as a python string.

  2. If the problem is known to be an LP problem, we can replace pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING with pywraplp.Solver.GLOP_LINEAR_PROGRAMMING or with pywraplp.Solver.CLP_LINEAR_PROGRAMMING. (Google Glop is an LP solver based on the simplex method developed by Google.)

  3. If the problem were to minimize, change solver.Maximize to solver.Minimize.

The following shows a programmatic way to solve \[\begin{array}{rl} \min & \mathbf{c}^\mathsf{T} \mathbf{x} \\ \mbox{s.t.} & \mathbf{A}\mathbf{x} = \mathbf{b} \\ & \mathbf{x} \geq \mathbf{0}. \end{array}\] 


from ortools.linear_solver import pywraplp

prob_name = 'another_example'
solver = pywraplp.Solver(prob_name,
                         pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

infinity = solver.infinity()
tol = solver.tolerance()

A = [[1., 1., 1.], [1., -1., 2]]
b = [2.,1.]
c = [-2.5, 3.4, 1.7]

m = len(A)
n = len(c)

x = [solver.NumVar(0., infinity, 'x' + str(i)) for i in range(n)]

for i in range(m):
    solver.Add( sum( [A[i][j] * x[j] for j in range(n)]) == b[i])
    
solver.Minimize( sum( [c[j]*x[j] for j in range(n)]))

result_status = solver.Solve()

if (result_status == pywraplp.Solver.OPTIMAL):
    print('Solution:')
    print('Objective value = ', solver.Objective().Value())
    for i in range(n):
        if x[i].solution_value() > 0.0:
            print( x[i].name() + ' = ', x[i].solution_value())
else:
    print("No optimal solution found.")

# The following exports the model as a string in LP format
model_as_string = solver.ExportModelAsLpFormat(False)
print( model_as_string )

# Save the string in a file
lp_file = open( prob_name + '.lp', 'w')
lp_file.write( model_as_string )
lp_file.close()

NEOS server for optimization

If one does not want to download and install any solvers for solving difficult mixed-integer linear programming problems (and many other types of optimization problems), one can use the NEOS server for optimization. There are many solvers to choose from. To solve a mixed-integer linear programming problem using IBM ILOG CPLEX, one can submit a file in CPLEX LP format here. The CPLEX LP format is described next.

CPLEX LP file format

A description of the CPLEX LP file format can be found here. However, it is perhaps easiest to illustrate with some examples. Consider \[\begin{array}{rrcrcrcrcl} \min & x_1 & - & 2x_2 & +& x_3 & & & \\ \mbox{s.t.} & x_1 & - & 2x_2 & + & x_3 & - & x_4 & \geq & 1 \\ & x_1 & + & x_2 & + & 2x_3 & - & 2x_4 & = & 3 \\ & & - & x_2 & + & 3 x_3 & - & 5x_4 & \leq & 7 \\ & & & & & x_3 & , & x_4 & \geq & 0. \end{array}\] The problem can be specified in LP format as follows: 


min x_1 - 2x_2 + x_3
st
 x_1 - 2x_2 +  x_3 -  x_4 >= 1
 x_1 +  x_2 + 2x_3 - 2x_4  = 3
     -  x_2 + 3x_3 - 5x_4 <= 7
bounds
 x_1 free
 x_2 free
end
Any variable that does not appear in the BOUNDS section is automatically assumed to be nonnegative. In other words, x_3 and x_4 above are automatically treated as bounded below by 0.

Mixed-integer linear programming problems can also be specified. Consider \[\begin{array}{rrcrcl} \max & x & - & 2y \\ \mbox{s.t.} & x & - & y & = & 2 \\ & 2x & + & 3y & \leq & 1 \\ & x & , & y & \geq & 0 \\ & & & y & \in & \mathbb{Z}. \end{array}\] The problem can be specified in LP format as follows: 


max x - 2y
st
  x -  y  = 2
 2x + 3y <= 1
general
 y
end

Note that we do not need a BOUNDS section because both variables are required to be nonnegative. However, \(y\) is an integer variable. Thus, it appears in the GENERAL section. In general, any variable that is required to be an integer must appear in the GENERAL section. If a variable is known to be binary, one can list it in the BINARY section instead of the GENERAL. For example, the problem \[\begin{array}{rrcrcl} \max & x & - & 2y \\ \mbox{s.t.} & x & - & y & = & 2 \\ & 2x & + & 3y & \leq & 1 \\ & & & y & \in & \{0,1\} \end{array}\] can be specified in two different ways: 


max x - 2y
st
  x -  y  =  2
 2x + 3y <=  1
bounds
 x free
 y <= 1
general
 y
end
 

max x - 2y
st
  x -  y  =  2
 2x + 3y <=  1
bounds
 x free
binary
 y
end

SoPlex

A special mention needs to be made to SoPlex, an open-source linear programming solver. It is free for academic use. Even though it does not support mixed-integer linear programming problems, it can return exact rational solutions whereas all the solvers mentioned above only return solutions as floating-point numbers. Binaries for Mac OS X and Windows are readily available for download.

Suppose that the problem for the example at the beginning is saved in CPLEX LP format in a plain-text file named eg.lp. The following is the output of running SoPlex in a Mac OS X command-line terminal: 

bash-3.2$ ./soplex-2.2.1.darwin.x86_64.gnu.opt -X --solvemode=2 -f=0 -o=0 eg.lp
SoPlex version 2.2.1 [mode: optimized] [precision: 8 byte] [rational: GMP 6.0.0] [githash: 267a44a]
Copyright (c) 1996-2016 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)

int:solvemode = 2
real:feastol = 0
real:opttol = 0

Reading (real) LP file  . . .
Reading took 0.00 seconds.

LP has 3 rows 4 columns and 11 nonzeros.

Initial floating-point solve . . .
Simplifier removed 0 rows, 0 columns, 0 nonzeros, 0 col bounds, 0 row bounds
Reduced LP has 3 rows 4 columns 11 nonzeros
Equilibrium scaling LP

type |   time |   iters | facts |  shift   |violation |    value
  L  |    0.0 |       0 |     1 | 4.00e+00 | 2.00e+00 | 0.00000000e+00
  E  |    0.0 |       1 |     2 | 0.00e+00 | 4.00e+00 | 3.00000000e+00
  E  |    0.0 |       2 |     3 | 0.00e+00 | 0.00e+00 | 1.00000000e+00

Floating-point optimal.
Max. bound violation = 0
Max. row violation = 1/4503599627370496
Max. reduced cost violation = 0
Max. dual violation = 0
Performing rational reconstruction . . .
Tolerances reached.
Solved to optimality.

SoPlex status       : problem is solved [optimal]
Solving time (sec)  : 0.00
Iterations          : 2
Objective value     : 1.00000000e+00


Primal solution (name, value):
x_1	              7/3
x_2	              2/3
All other variables are zero.
bash-3.2$

The option -X asks the solver to display the primal rational solution. The options --solvemode=2 invokes iterative refinement for solving for a rational solution. The options -f=0 -o=0 set the primal feasibility and dual feasibility tolerances to 0. Without these options, one might get only approximate solutions to the problem. If we remove the last three options and replace -X with -x, we obtain the following instead: 

bash-3.2$ ./soplex-2.2.1.darwin.x86_64.gnu.opt -x eg.lp
SoPlex version 2.2.1 [mode: optimized] [precision: 8 byte] [rational: GMP 6.0.0] [githash: 267a44a]
Copyright (c) 1996-2016 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)

Reading (real) LP file  . . .
Reading took 0.00 seconds.

LP has 3 rows 4 columns and 11 nonzeros.

Simplifier removed 0 rows, 0 columns, 0 nonzeros, 0 col bounds, 0 row bounds
Reduced LP has 3 rows 4 columns 11 nonzeros
Equilibrium scaling LP
type |   time |   iters | facts |  shift   |violation |    value
  L  |    0.0 |       0 |     1 | 4.00e+00 | 2.00e+00 | 0.00000000e+00
  E  |    0.0 |       1 |     2 | 0.00e+00 | 4.00e+00 | 3.00000000e+00
  E  |    0.0 |       2 |     3 | 0.00e+00 | 0.00e+00 | 1.00000000e+00
 --- transforming basis into original space
  L  |    0.0 |       0 |     1 | 0.00e+00 | 0.00e+00 | 1.00000000e+00
  L  |    0.0 |       0 |     1 | 0.00e+00 | 0.00e+00 | 1.00000000e+00

SoPlex status       : problem is solved [optimal]
Solving time (sec)  : 0.00
Iterations          : 2
Objective value     : 1.00000000e+00


Primal solution (name, value):
x_1	      2.333333333
x_2	      0.666666667
All other variables are zero (within 1.0e-16).
bash-3.2$

There are many solver options that one can specify. To get the list of all the options, simply run the solver without options and arguments.

Worked examples

  1. Use the IBM ILOG CPLEX solver at the NEOS server to solve the following: \begin{align*} \min \ & x + y + z\\ \text{s.t.} \ & 7x - 2y - z \geq 5 \\ & -x + 3y + 4z \geq 1 \\ & x, y, z \geq 0. \end{align*}

  2. Use Wolfram Alpha to solve \begin{align*} \min \ & x + 2y + 3z\\ \text{s.t.} \ & 3x - 2y + z \geq 2 \\ & -2x + y + 5z \geq 3 \\ & x, y, z \geq 0. \end{align*}