Thursday, July 29, 2010

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.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);
// 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:
2. 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.
3. We must use the primal simplex algorithm. If we comment out that line (and restore the previously removed lines), we get the following output:
4. 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.

5. 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.
6. 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).

1. This is a good topic and often misunderstood. I never realized why cplex handles infeasibility/unboundeness in this way. In my opinion this is the way to do it :

where a farkas certificate is always returned. You can check the certificate and how "strong" it is.

2. I should say I used to work for mosek, so I might be a bit biased (though I am convinced this is the right way ;-)).

3. @Bo: Thanks for the reminder about Farkas certificates. CPLEX does provide a Farkas certificate via IloCplex.dualFarkas, but I think only for infeasible primals.

4. Hi Paul,

In CPLEX, is it possible to find more than one extreme ray from a single solve? I have a model where there are many infeasibilities. It seems CPLEX provides one extreme ray at a time per solve. If I fix one infeasibility then I could resolve to get another one.

Regards,
Vivek.

1. Vivek: You cannot get additional rays just from a function call. CPLEX stops as soon as it recognizes that the problem is unbounded, and it only knows the ray that allowed it to recognize the unboundedness.

I just wrote a new post describing one way to find some (but not all) additional rays. It requires modifying the objective function of the unbounded problem and re-solving.

Cheers,
Paul

2. Thanks Paul!

5. could you provide us with mathematical formula how to compute the extreme direction, and do you have a good resource of linear programming (contains a lot of problem set )

1. I don't think there is a closed form formula (equation) for an extreme direction. A recession direction for $Ax \ge b, \, x \ge 0$ will be any direction $d \ge 0$ for which $Ad \ge 0$. Which recession directions are extreme is a bit harder to characterize algebraically.

The usual approach to finding an extreme direction is to apply one of the simplex algorithms to the LP. If it has recession directions along which the objective improves, at some point the algorithm will terminate with a message saying something equivalent to "I know where I want to pivot, but there's no limit to how far I can go on the new edge". The specific condition depends on which simplex algorithm (primal, dual, primal-dual, ...) you are using. Note that recession directions along which the objective function gets worse are generally not detected (and generally of no interest).

Regarding references, I've seen (and used) books I thought were rather good, with adequate end-of-chapter exercises (but no solutions). They're mostly out of print now, since I last taught LP many, many years ago. You might ask for recommendations from younger people on OR-Exchange (link in the right margin).

If this is your first time commenting on the blog, please read the Ground Rules for Comments. In particular, if you want to ask an operations research-related question not relevant to this post, consider asking it on OR-Exchange.