I’m teaching undergrad linear optimization this semester.  Mathematica can be helpful for visualizing LPs.  It also has some LP solving capabilities.  (Special thanks on this post goes to my colleague Russ Rushmeier, with whom I’m teaching the course, for the example LP, some of the commands, and great discussions.)

Here is the example LP we’ll use:

max z = x1 + x2 + x3

subject to:

3 x1 + 2 x2 + x3 ≤ 6
x1 – x2 + x3 ≤ 1
x1 + 2 x2 +  x3 ≤ 4
x1 ≥ 0, x2 ≥ 0, x3 ≥ 0.

First we’ll plot the feasible region, which is 3-dimensional, using the “RegionPlot3D” command:

fr = RegionPlot3D[3 x1 + 2 x2 + x3 <= 6 &&
x1 – x2 + x3 <= 1 &&
x1 + 2 x2 + 1 x3 <= 4 &&
x1 >= 0 && x2 >= 0 && x3 >= 0,
{x1, 0, 3}, {x2, 0, 3}, {x3, 0, 3},
Mesh -> All, PlotPoints -> 100, AxesLabel -> Automatic]

In Mathematica you can rotate it around to get a feel for the polyhedral shape.  For 3-variable problems, constraints are represented as planes dividing 3-space into two parts.

We can visualize the objective contours as planes.  Using the objective z = x1 + x2 + x3, we can plot the contours z = 1.5, z = 2.5 (solving each equation for x3 and plotting the surface):

zcontours = Plot3D[{1.5 – x1 – x2, 2.5 – x1 – x2} , {x1, 0, 3}, {x2, 0, 3},AxesLabel -> Automatic];
Show[fr, zcontours]

(The feasible region and z contours are each stored in separate graphics objects that are then displayed together using the “Show” command.)

The contour method of solving an LP is usually illustrated with a 2 variable problem, but with these graphics set up in Mathematica, it is not hard to see it in 3D.  Instead of displaying 2 static contours, we keep one fixed and allow the other to vary by using the “Manipulate” command or the “Animate” command.  Animate goes by itself and is what I used to generate this animated gif.  Manipulate lets you control the motion with a slider.  Here I’ve allowed the second contour to move between z values of 1 and 5 and you can see (hopefully) it pass through optimal solutions along the way.  (The 2 extreme point optimal solutions are indicated with blue dots.)

Manipulate[Show[fr, Plot3D[{1.5 – x1 – x2, z – x1 – x2} , {x1, 0, 3}, {x2, 0,3}, AxesLabel -> Automatic]], {z, 1, 5}]

To solve the LP, we first try Mathematicas “Maximize” command.  The syntax is pretty straightforward – objective first followed by &&-separated constraints including sign restrictions, then ending with a list of the variables to be found:

Maximize[x1 + x2 + x3, 3 x1 + 2 x2 + x3 <= 6 &&
x1 –    x2 + x3 <= 1  &&
x1 + 2 x2 + x3 <= 4 &&
x1 >= 0 &&  x2 >= 0 && x3 >= 0, {x1, x2, x3}]

Note it reports the lower of the two optimal extreme points but no indication that alternate optimal solutions exist.  Reading the help file for Maximize shows many other capabilities, but I did not notice anything addressing this issue.

A better option may be the “LinearProgramming” command.  This assumes a min objective (the first list), followed by the coefficients of the A matrix, and then the values of the b vector where the -1 next to each value serves as a flag indicating a less-than-or-equal-to constraint:

LinearProgramming[{-1, -1, -1}, {{3, 2, 1}, {1, -1, 1}, {1, 2, 1}}, {{6, -1}, {1, -1}, {4, -1}}]

The same solution is found, but again we can’t tell if it’s the only one.

LinearProgramming appears to use the Simplex method by default but can use Interior Point methods if specified:

LinearProgramming[{-1, -1, -1}, {{3, 2, 1}, {1, -1, 1}, {1, 2, 1}}, {{6, -1}, {1, -1}, {4, -1}}, Method -> “InteriorPoint”]

So now we see evidence of alternate optimals.  This solutions is on the interior of the optimal edge, so a little more work would be required to fully specify all optimal solutions.

We handle this differently in the course.  At this point, we are learning about how the simplex method works and seeing the (possible) existence of alternate optimal solutions through zero objective row coefficients in the simplex tableau.  From there we pivot to an adjacent optimal extreme point, repeating if necessary.  Then we can express the set accordingly.  Soon we will move into other technology to see how this plays out, starting with Excel Solver.  But the visualization capabilities of Mathematica for 3D LPs do help illustrate what is going on.