Linear programming in Python

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:

import pymathprog as pm
pm.begin('furniture')
# pm.verbose(True)
x, y = pm.var('x, y') # variables
pm.maximize(30 * x + 60 * y, 'profit')
0.5*x + 1.5*y <= 30 # wood
15*x + 20*y <= 600 # laminate
pm.solve()
print('\nMax profit:',pm.vobj())
pm.sensitivity()
pm.end()

with output:

GLPK Simplex Optimizer, v4.65
2 rows, 2 columns, 4 non-zeros
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
*     2: obj =   1.440000000e+03 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND

Max profit: 1440.0

PyMathProg 1.0 Sensitivity Report Created: 2018/10/28 Sun 21:42PM
================================================================================
Variable            Activity   Dual.Value     Obj.Coef   Range.From   Range.Till
--------------------------------------------------------------------------------
*x                        24            0           30           20           45
*y                        12            0           60           40           90
================================================================================
================================================================================
Constraint       Activity Dual.Value  Lower.Bnd  Upper.Bnd RangeLower RangeUpper
--------------------------------------------------------------------------------
 R1                    30         24       -inf         30         20         45
 R2                   600        1.2       -inf        600        400        900
================================================================================

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:

from scipy.optimize import linprog
linprog([-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

    fun: -1440.0
message: 'Optimization terminated successfully.'
    nit: 2
  slack: array([0., 0.])
 status: 0
success: True
      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.

 
comments powered by Disqus