I’ve been away helping Bain set up a new capability in generative AI, on the back of our OpenAI partnership, and am coming back to standard data science problems. In this post, I’d like to share some experiments doing linear programming optimization, using Gonum’s optimization library. It’s a very basic facility, but has helped me learn how to set up problems using the venerable Simplex algorithm, first developed by George Danzig in 1947.

First of all, I wouldn’t recommend using this approach for real problem solving, since it’s rather cumbersome, as you will see. Instead, use Pulp or another Python library, or a dedicated mathematical programming tool such as GLPK (open source) or Gurobi (commercial license). Nevertheless, I found this approach helpful for finally understanding how problems are represented as matrices, and it serves as a readily available, minimalist approach to solving simpler linear optimization problems.

A simple problem

Here is a simple example LP problem, from Brunel University:

A company is involved in the production of two items (X and Y). The resources need to produce X and Y are twofold, namely machine time for automatic processing and craftsman time for hand finishing. The table below gives the number of minutes required for each item:

         Machine time Craftsman time
Item X   13           20
     Y   19           29

The company has 40 hours of machine time available in the next working week but only 35 hours of craftsman time. Machine time is costed at £10 per hour worked and craftsman time is costed at £2 per hour worked. Both machine and craftsman idle times incur no costs. The revenue received for each item produced (all production is sold) is £20 for X and £30 for Y. The company has a specific contract to produce 10 items of X per week for a particular customer.

Formulating the problem

A simple formulation to this problem into decision variables, objective function, and constraints, looks as follows:

Let
    x be the number of items of X
    y be the number of items of Y

Maximise
    20x + 30y - 10(machine time worked) - 2(craftsman time worked)

Subject to:
    13x + 19y <= 2400       // machine time
    20x + 29y <= 2100       // craftsman time
    x >= 10                 // contract
    x,y >= 0

The objective function can be expanded and simplified as:

Maximise
    20x + 30y - 10(13x + 19y)/60 - 2(20x + 29y)/60
    => 17.1667x + 25.8667y

Converting to matrix

The above formulation can be input almost verbatim into the various standard optimization libraries (example below). To use Gonum’s simplex optimizer, however, you need to convert it to a matrix. This process is not described in the documentation, and I could not find any pages that show examples with explanation. So here we go.

The call to the optimizer will look like this:

opt, x, err := lp.Simplex(c, A, b, 0, nil)

There are three important arguments to Simplex:

• c is a vector (array of float64 numbers) of the objective function coefficients (negative to minimize)
• A is a matrix of the left hand side coefficients for the constraints, one row per constraint
• b is a vector of the right hand side values of the constraints

Let’s take them in turn.

Objective function coefficients

The objective function is

Maximise:
    17.1667x + 25.8667y

There are only two decision variables (x and y), so we are only interested in two columns here. However, we need to add three more columns for the “slack” variables, one per constraint, otherwise there will be no feasible solution if the optimum values do not exactly equal the constraints.

Note: can we make do with two slack variables, one per decision variable (x and y) instead of one for each of the three constraints?

Furthermore, the Simplex algorithm does minimization, so to maximize, we need to multiply the coefficents by -1.

So c, our vector for the objective function coefficients, is:

c := []float64{-17.2, -25.9, 0, 0, 0}

LHS constraint coefficients

A is a Gonum matrix of the left hand side coefficients for the constraints. There is one row per constraint. There is also one column per decision variable (x and y), but also additional columns with diagonal ones, one per constraint.

Note that

• <= constraints are converted to = by adding slack variables
• <= constraints are represented with the same sign as the problem definition
• Any >= constraints need to have the coefficients multiplied by -1, both in the left hand side (matrix A), and also the vector of right-hand side values (c, see below)

Our constraints are:

Subject to:
    13x + 19y <= 2400       // machine time
    20x + 29y <= 2100       // craftsman time
    x >= 10                 // contract
    x,y >= 0

The matrix is 3 rows x 5 columns (note that x,y >= 0 does not need to expressed, as the algorithm takes care of this constraint):

13 19  1  0  0
20 29  0  1  0
-1  0  0  0  1

The first row is 13x + 19y, the second row is 20x + 20y, and the third row is 1x (converted to negative because it is >= instead of <=).

The Go matrix (using Gonum’s mat module) is:

A := mat.NewDense(3, 5, []float64{
    13, 19, 1, 0, 0,
    20, 29, 0, 1, 0,
    -1, 0, 0, 0, 1})

Right hand constraint values

Finally, b is a vector of the right hand side values of the constraints, the bounds to which we are restricting the various constraints.

Again, given these constraints:

Subject to:
    13x + 19y <= 2400       // machine time
    20x + 29y <= 2100       // craftsman time
    x >= 10                 // contract

The bounds are 2400, 2100, and 10. But the 10 needs to be -10, since it is a >= instead of <= constraint. So the Go array is:

b := []float64{2400, 2100, -10}

Putting it all together

The complete program looks as follows:

package main

import (
    "fmt"
    "gonum.org/v1/gonum/mat"
    "gonum.org/v1/gonum/optimize/convex/lp"
)

func main() {

    // Problem formulation
    c := []float64{-17.2, -25.9, 0, 0, 0}
    A := mat.NewDense(3, 5, []float64{
        13, 19, 1, 0, 0,
        20, 29, 0, 1, 0,
        -1, 0, 0, 0, 1})
    b := []float64{2400, 2100, -10}

    // Run Simplex algorithm, last parameter is initialBasic
    opt, x, err := lp.Simplex(c, A, b, 0, nil)
    if err != nil {
        fmt.Println(err.Error())
        return
    }

    fmt.Printf("Optimum = %f, x = %f, y = %f\n", opt, x[0], x[1])
}

Note that in addition to c, A, and b, lp.Simplex also requires a number as the tolerance (0 as in the example above), and optionally a vector representing the starting values.

Also note that lp.Simplex returns the optimum value, an array of coefficients, and an error (e.g., solution is infeasible, or problem is unbounded). The array of solution coefficients corresponds to the columns in the A matrix, so you only need to the first two in the example above (altough the others are interesting for seeing the slack values).

You will need to install Gonum’s matrix and simplex modules:

go get gonum.org/v1/gonum

So to build and run the program, do the following:

mkdir lpexample
cd lpexample
vim main.go
(paste the above code and save)
go mod init lpexample
go get gonum.org/v1/gonum
go build
./lpexample

You should get the solution x = 10, y = 65.52, max = 1866.5

Formulation in Python

As mentionned above, such programs are much easier to create and understand using a Python library or a dedicated mathematical programming tool. For comparison, below is the same problem formulated in Python using Pulp. It’s a bit longer, but easier to read and maintain, because you don’t need to keep track of the columns and signs in the matrices and vectors. It also allows you to switch between different solvers.

from pulp import *

# Create the 'prob' variable to contain the problem data
prob = LpProblem("Example", LpMaximize)

# Problem variables (use LpInteger for integer solution)
x = LpVariable("x", 0, None, LpContinuous)
y = LpVariable("y", 0, None, LpContinuous)

# Objective function
prob += 17.2*x + 25.9*y

# Constraints
prob += 13*x + 19*y <= 2400 # machine time constraint
prob += 20*x + 29*y <= 2100 # labor time constraint
prob += x >= 10 # contract commitment

# Solve the problem
status = prob.solve()

# Print results
print("Status:", LpStatus[status])
print("x:", value(x))
print("y:", value(y))
print("Objective:", value(prob.objective))

Again, you should get the solution x = 10, y = 65.52, max = 1866.5 (or slightly different if you chose to use LpInteger for integer programming).