Finding Extreme Rays in CPLEX

It's funny how questions in technical forums sometimes cluster. After going quite some time without seeing any mention of extreme rays on any of the CPLEX-related forums (nor on sci.op-research, nor OR-Exchange), in the past month or so I've encountered at least three threads relating to the topic. So perhaps it's fodder for a blog post.

Some terminology first. If a convex polyhedron is unbounded, it has one or more recession directions, directions of rays that remain within the polyhedron. The set of recession directions is convex, so “one or more” means “one or uncountably many”. Nifty properties include that a recession direction from one point in the polyhedron is a recession direction from every point in the polyhedron, and that the recession directions form a convex cone. An extreme direction is a recession direction that cannot be written as a non-trivial convex combination of other recession directions; the corresponding extreme ray is a boundary of the cone of feasible rays.

Now let's move to the realm of linear programs. If a linear program is unbounded, then the feasible region (a convex polyhedron) must be unbounded in a direction that improves the objective function. The primal simplex method recognizes this when it reaches a vertex at which there is a variable with favorable reduced cost to pivot out of the basis (indicating an edge on which to move) but there is no candidate to enter the basis (because no constraint becomes binding as you move along the edge). The edge in question is an extreme ray anchored at an extreme point of the feasible region.

Consider the problem\begin{eqnarray*} \mathrm{maximize} & 3x+2y\\ \mathrm{s.t.} & x-2y & \le1\\ & -2x+y & \le1\\ & x+y & \ge2\\ & x,y & \ge0\end{eqnarray*}whose feasible region is shown in Figure 1 below.

There are two extreme rays, one anchored at $A=\left(\frac{5}{3},\frac{1}{3}\right)$ with direction $\left(2,1\right)$ and the other anchored at $B=\left(\frac{1}{3},\frac{5}{3}\right)$ with direction $\left(1,2\right)$. I'll code it in CPLEX using the Java API (other APIs will behave similarly), omitting imports, try/catch blocks etc.

IloCplex cplex = new IloCplex();
IloNumVar x = cplex.numVar(0, Double.MAX_VALUE, "x");
IloNumVar y = cplex.numVar(0, Double.MAX_VALUE, "y");
cplex.addMaximize(cplex.sum(cplex.prod(3.0, x),
                  cplex.prod(2.0, y)));
cplex.addLe(cplex.sum(x, cplex.prod(-2.0, y)), 1.0);
cplex.addLe(cplex.sum(y, cplex.prod(-2.0, x)), 1.0);
cplex.addGe(cplex.sum(x, y), 2.0);
// turn off the presolver
cplex.setParam(IloCplex.BooleanParam.PreInd, false);
// select the primal simplex method
cplex.setParam(IloCplex.IntParam.RootAlg,
               IloCplex.Algorithm.Primal);
// solve the model
cplex.solve()
// verify primal is unbounded
System.out.println("CPLEX status: "
                   + cplex.getCplexStatus());
// get the current solution
System.out.println("Last extreme point: x = "
                   + cplex.getValue(x)
                   + ", y = "
                   + cplex.getValue(y));
// get a ray
IloLinearNumExpr ray = cplex.getRay();
System.out.println("getRay returned " + ray.toString());
// parse it
IloLinearNumExprIterator it = ray.linearIterator();
double xcoef = 0, ycoef = 0;
while (it.hasNext()) {
  IloNumVar v = it.nextNumVar();
  if (v.equals(x)) {
    xcoef = it.getValue();
  }
  else if (v.equals(y)) {
    ycoef = it.getValue();
  }
}
System.out.println("Extreme ray direction = ("
                   + xcoef + ", " + ycoef + ")");

The output is as follows:

The output is as follows:
Iteration log . . .

Iteration:     1    Infeasibility =             1.000000

Switched to steepest-edge.

CPLEX status: Unbounded

Last extreme point: x= 1.6666666666666665,
                    y = 0.3333333333333335

getRay returned 0.6666666666666666*x
                + 0.33333333333333337*y

Extreme ray direction = (0.6666666666666666,
                         0.33333333333333337)

Note that the primal simplex method stopped at vertex $A$ and the direction vector is within a scalar multiple of $(2,1)$, so it is the correct direction for the extreme ray there.

Things to note:

1. We must turn off presolving to find the ray. If we comment out that line, the attempt to access the vertex causes an exception:

CPLEX Error  1217: No solution exists.

If we also comment out the lines that try to get the vertex values (seeking just the ray and not its anchor point), the call to getRay() throws an exception:

CPLEX Error  1254: Unbounded solution required.

2. We must use the primal simplex algorithm. If we comment out that line (and restore the previously removed lines), we get the following output:

Last extreme point: x = 0.0, y = 0.0
getRay returned 1.0*x + 1.0*y
Extreme ray direction = (1.0, 1.0)

Note that the anchor point $(0,0)$ is not actually a vertex of the feasible region (and in fact is not in the feasible region), and the direction vector $(1,1)$, while a legitimate recession direction, is not an extreme recession direction.


3. Our arrival at vertex $A$ rather than vertex $B$ was due to the choice of the objective function. Tilt the objective function so that its normal vector is a bit closer to vertical, and we end up at $B$ instead.
4. The bad news is that in most APIs we need to fiddle with an iterator to parse the direction vector, since getRay() returns a linear expression rather than a double[]API, CPXgetray() passes the ray coefficients in a double[] argument, including any zeros (no attempt to exploit any sparsity).