Few days ago I decided to create a C++ implementation of the Simplex algorithm, which is a tool for solving Linear Programming problems. One of my motivations could have been my recent encounter with the Eigen linear algebra library, which really surprised me with its beautiful syntax and easy usage. So I wanted to gain a little experience with Eigen, although some factorization-related project might have suited this purpose better. Another motivational factor was a not yet published project of mine, which might get finished in 1 or 2 months, and would give LP a real application.
Although in theory the Simplex method is quite simple, but its implementation is much harder and “dirtier” than one would think. I coded only its simplest form – which uses dense matrices – and yet found myself in the middle of complications and rare cases. None of them was mentioned in the books and articles I read. Originally I wanted to undertake its sparse matrix version, which maintains an LU decomposition, but after reading the complaining comments of developers having 10-20 years experience on this field, I’ve changed my mind quickly. Still this little trip allowed me a deeper insight into the concept of linear programming.
The source files are available on GitHub:
Two example problems are solved in the main.cpp above. Both are 2 dimensional, although the class can work with arbitrary matrix sizes, which fit into memory. One of the implemented problems is maximization while the other is minimization.
The maximization problem can be formulated in matrix form like:
Where x contains variables, which with the coefficients of A and b forms a set of linear inequalities, called: constraints. The c matrix contains the coefficients of the objective function cx, which we want to maximize. This is only a 2 dimensional problem, so a simple graphical representation is possible: (Some renaming was done: x1 = x and x2 = y)
The blue area is defined by the linear inequalities, the task is to find the topmost point, where the red line touches the blue area, if we start to move it downwards, without changing its slope. Easy to see that the solution is the point: (5, 8)
Now its time for the minimization problem:
And the geometric representation is: (Renamed again: x1 = x and x2 = y)
See the worked out example code: https://github.com/bolner/SimplexSolver/blob/master/main.cpp
Notice that the A and b matrices are not passed separately to the solver class (which is usual in most simplex implementations) but they form one big constraint matrix, where b is the rightmost column, like in the algebraic inequality form.