Linear programming in Python

Oct 28, 2018 linear-programming python
Share on:

For my elementary linear programming subject, the students (who are all pre-service teachers) use Excel and its Solver as the computational tool of choice. We do this for several reasons: Excel is software with which they're likely to have had some experience, also it's used in schools; it also means we don't have to spend time and mental energy getting to grips with new and unfamiliar software. And indeed the mandated curriculum includes computer exploration, using either Excel Solver, or the Wolfram Alpha Linear Programming widget.

This is all very well, but I balk at the reliance on commercial software, no matter how widely used it may be. And for my own exploration I've been looking for an open-source equivalent.

In fact there are plenty of linear programming tools and libraries; two of the most popular open-source ones are:

  • The GNU Linear Programming Kit, GLPK
  • Coin-or Linear Programming, Clp

There's a huge list on wikipedia which includes open-source and proprietary software.

For pretty much any language you care to name, somebody has taken either GLPK or Clp (or both) and produced a language API for it. For Python there's PuLP; for Julia there's JuMP; for Octave there's the `glpk` command, and so on. Most of the API's include methods of calling other solvers, if you have them available.

However not all of these are well documented, and in particular some of them don't allow sensitivity analysis: computing shadow prices, or ranges of the objective coefficients. I discovered that JuMP doesn't yet support this - although to be fair sensitivity analysis does depend on the problem being solved, and the solver being used.

Being a Python aficionado, I thought I'd check out some Python packages, of which a list is given at an operations research page.

However, I then discovered the Python package PyMathProg which for my purposes is perfect - it just calls GLPK, but in a nicely "pythonic" manner, and the design of the package suits me very well.

A simple example

Here's a tiny two-dimensional problem I gave to my students:

A furniture workshop produces chairs and tables. Each day 30m2 of wood board is delivered to the workshop, of which chairs require 0.5m2 and tables 1.5m2. (We assume, of course, that all wood is used with no wastage.) All furniture needs to be laminated; there is only one machine available for 10 hours per day, and chairs take 15 minutes each, tables 20 minutes. If chairs are sold for $30 and tables for $60, then maximize the daily profit (assuming that all are sold).

Letting \(x\) be the number of chairs, and \(y\) be the number of tables, the problem is to maximize \[ 30x+60y \] given

\begin{align*} 0.5x+1.5y&\le 30\\ 15x+20y&\le 600\\ x,y&\ge 0 \end{align*}

Problems don't get much simpler than this. In pyMathProg:

 1import pymathprog as pm
 2pm.begin('furniture')
 3# pm.verbose(True)
 4x, y = pm.var('x, y') # variables
 5pm.maximize(30 * x + 60 * y, 'profit')
 60.5*x + 1.5*y <= 30 # wood
 715*x + 20*y <= 600 # laminate
 8pm.solve()
 9print('\nMax profit:',pm.vobj())
10pm.sensitivity()
11pm.end()

with output:

 1GLPK Simplex Optimizer, v4.65
 22 rows, 2 columns, 4 non-zeros
 3*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
 4*     2: obj =   1.440000000e+03 inf =   0.000e+00 (0)
 5OPTIMAL LP SOLUTION FOUND
 6
 7Max profit: 1440.0
 8
 9PyMathProg 1.0 Sensitivity Report Created: 2018/10/28 Sun 21:42PM
10================================================================================
11Variable            Activity   Dual.Value     Obj.Coef   Range.From   Range.Till
12--------------------------------------------------------------------------------
13*x                        24            0           30           20           45
14*y                        12            0           60           40           90
15================================================================================
16================================================================================
17Constraint       Activity Dual.Value  Lower.Bnd  Upper.Bnd RangeLower RangeUpper
18--------------------------------------------------------------------------------
19 R1                    30         24       -inf         30         20         45
20 R2                   600        1.2       -inf        600        400        900
21================================================================================

From that output, we see that the required maximum is $1440, obtained by making 24 chairs and 12 tables. We also see that the shadow prices for the constraints are 24 and 1.2. Furthermore, the ranges of objective coefficients which will not affect the results are \([20,45]\) for prices for chairs, and \([40,90]\) for table prices.

This is the simplest API I've found so far which provides that sensitivity analysis.

Note that if we just want a solution, we can use the linprog command from scipy:

1from scipy.optimize import linprog
2linprog([-30,-60],A_ub=[[0.5,1.5],[15,20]],b_ub=[30,600])

linprog automatically minimizes a function, so to maximize we use a negative function. The output is

1    fun: -1440.0
2message: 'Optimization terminated successfully.'
3    nit: 2
4  slack: array([0., 0.])
5 status: 0
6success: True
7      x: array([24., 12.])

The negative value given as fun above simply reflects that we are entering a negative function. In respect of our problem, we simply negate that value to obtain the required maximum of 1440.