Advanced MIP Start
CPLEX allows you to use a heuristic to generate an initial incumbent integer solution. Recall that at each node (for minimization problems), if the LP relaxation takes a value greater than the current incumbent, then we can prune the node without further branching. See the figure below:
Thus having a good incumbent initial incumbent solution will reduce the total number of nodes you visit in branch and bound. Assuming you can generate your incumbent quickly, this could reduce the running time.
Advanced MIP Start in CPLEX
The class IloCplex provides a method for setting the advanced start:
Method Name |
Return Type |
Arguments |
Description |
|---|---|---|---|
addMIPStart |
int |
IloNumVar[] vars, double[] values |
Adds the MIP start specified by its variables and values to the current problem. Multiple calls to this function will overwrite previous calls, not add additional starting points. Return indicates?? |
The advanced MIP start actually has a lot of features that we aren't using, as explained in the documentation. It can accept only partial or infeasible solutions and then will try and "repair" them to get an incumbent. It can also apply local search techniques to improve the quality of the incumbent, but we will be doing that with a custom algorithm soon instead.
Christofides Algorithm
Christofides Algorithm is an approximation algorithm that gives a 1.5 approximation factor for metric TSP problems (TSP problems where distance/cost function between nodes is a a metric). In fact, the algorithm can be run on any TSP problem on a complete graph, although we have performance guarantee if the distances are not a metric. As all of the problems we are considering use the Euclidean metric for distance (see also Euclidean TSP), Christofides algorithm should give us good performance.
Implementing Advanced MIP Start for TspIpSolver
We have provided an implementation of Chirstofides algorithm for your convenience in the class ChristofidesSolver. The solver has a single method:
Method Name |
Return Type |
Arguments |
Description |
|---|---|---|---|
approximateBestTour |
List<E> tour |
Set<E> suggestedEdges |
Returns the tour given by running Christofides algorithm on the input graph and edge weights. The algorithm attempts to include all edges in suggestedEdges by adjusting the weights when forming a MST (if suggested edges is empty, the result is classical Christofides). If suggestedEdges contains a loop, will return null instead of a tour. |
To incorporate the solver, add the field to TspIpSolver
private ChristofidesSolver<V,E> christofides;
We need a utility method that given a set of edges, will create an array of zeros and ones. Directed after edgesUsed() in TspIpSolver, and the method
/**
* Given a subset of the edges, produces a binary vector indicating which edges were included using the same ordering as edgeVariablesAsArray
* @param edgesUsed a subset of the set of edges in the grpah
* @return an array of zeros and ones where the ith entry corresponds to the ith edge variable from edgeVariablesAsArray
*/
private double[] inverseEdgesUsed(Set<E> edgesUsed){
double[] edgeVals = new double[this.edgeVariablesAsArray.length];
for(int i =0; i < edgeVals.length; i++){
edgeVals[i] = edgesUsed.contains(edgeVariables.inverse().get(edgeVariablesAsArray[i])) ? 1 : 0;
}
return edgeVals;
}
Finally, at the very end of the constructor, add the lines
if(options.contains(Option.christofidesApprox)){
System.out.println("Beginning Christofides...");
ChristofidesSolver<V,E> christofides = new ChristofidesSolver<V,E>(graph,tspInstance.getEdgeWeights());
List<E> christofidesApproxTour = christofides.approximateBestTour(new HashSet<E>());
Set<E> submittedTour = new HashSet<E>(christofidesApproxTour);
double[] edgeVals = inverseEdgesUsed(submittedTour);
cplex.addMIPStart(edgeVariablesAsArray, edgeVals);
}
Profiling on TSPLIB
Again running d493, we see that Christofides gives improvements across the board
We also see that there is still considerable room for improvement, as the gap between Christofides and optimum is quite large (although much less than 3/2).
In the final step of Christofides algorithm, where a Eulerian cycle is converted to a Hamiltonian cycle, we make no attempt to optimize this process. While it was shown that optimal short cutting is NP hard, there are some greedy heuristics that give a performance boost (see here). Ultimately, we will take a different approach, and instead apply local search to the output of Christofides.
Complete Source
If all else fails, at this point, your code should look like this:
Toggle Complete Source
public class TspIpSolver<V,E> {
public static enum Option{
lazy,userCut, randomizedUserCut, christofidesApprox, christofidesHeuristic,twoOpt,incumbent;
}
private EnumSet<Option> options;
private IloCplex cplex;
private TspInstance<V,E> tspInstance;
private final ImmutableBiMap<E,IloIntVar> edgeVariables;
private ImmutableSet<E> edgesInOpt;
private double optVal;
private IloIntVar[] edgeVariablesAsArray;
private ChristofidesSolver<V,E> christofides;
public TspIpSolver(TspInstance<V,E> tspInstance) throws IloException{
this(tspInstance,EnumSet.of(Option.lazy, Option.userCut, Option.christofidesApprox, Option.christofidesHeuristic));
}
public TspIpSolver(TspInstance<V,E> tspInstance, EnumSet<Option> options) throws IloException{
this.options = options;
this.tspInstance = tspInstance;
this.cplex = new IloCplex();
UndirectedGraph<V,E> graph = tspInstance.getGraph(); //for convenience, we will be using this a lot
this.edgeVariables = Util.makeBinaryVariables(cplex, graph.getEdges());
edgeVariablesAsArray = edgeVariables.inverse().keySet().toArray(new IloIntVar[]{});
//the degree constraints
for(V vertex: graph.getVertices()){
cplex.addEq(Util.integerSum(cplex, edgeVariables,
graph.getIncidentEdges(vertex)), 2);
}
//the objective
cplex.addMinimize(Util.integerSum(
cplex, edgeVariables, graph.getEdges(),tspInstance.getEdgeWeights()));
if(options.contains(Option.lazy)){
cplex.use(new IntegerCutCallback());
}
if(options.contains(Option.userCut)){
cplex.use(new MinCutCallback());
if(options.contains(Option.randomizedUserCut)){
System.err.println("Warning, userCut and randomizedUserCut both selected, ignoring randomizedUserCut");
}
}
else if(options.contains(Option.randomizedUserCut)){
cplex.use(new RandomMinCutCallback());
}
if(options.contains(Option.christofidesApprox)){
System.out.println("Beginning Christofides...");
ChristofidesSolver<V,E> christofides = new ChristofidesSolver<V,E>(graph,tspInstance.getEdgeWeights());
List<E> christofidesApproxTour = christofides.approximateBestTour(new HashSet<E>());
Set<E> submittedTour = new HashSet<E>(christofidesApproxTour);
double[] edgeVals = inverseEdgesUsed(submittedTour);
cplex.addMIPStart(edgeVariablesAsArray, edgeVals);
}
}
public void solve() throws IloException{
if(!cplex.solve()){
throw new RuntimeException();
}
optVal = cplex.getObjValue();
edgesInOpt = ImmutableSet.copyOf(edgesUsed(cplex.getValues(edgeVariablesAsArray)));
cplex.end();
}
public ImmutableSet<E> getEdgesInOpt(){
return edgesInOpt;
}
public double getOptVal(){
return optVal;
}
/**
* assumes edgeVarVals[i] in {0,1}, up to some tolerance.
* @param edgeVarVals
* @return the edges where the variable takes the value one.
*/
private Set<E> edgesUsed(double[] edgeVarVals){
Set<E> ans = new HashSet<E>();
for(int e = 0; e < edgeVarVals.length; e++){
if(Util.doubleToBoolean(edgeVarVals[e])){
ans.add(edgeVariables.inverse().get(edgeVariablesAsArray[e]));
}
}
return ans;
}
/**
* Given a subset of the edges, produces a binary vector indicating which edges were included using the same ordering as edgeVariablesAsArray
* @param edgesUsed a subset of the set of edges in the grpah
* @return an array of zeros and ones where the ith entry corresponds to the ith edge variable from edgeVariablesAsArray
*/
private double[] inverseEdgesUsed(Set<E> edgesUsed){
double[] edgeVals = new double[this.edgeVariablesAsArray.length];
for(int i =0; i < edgeVals.length; i++){
edgeVals[i] = edgesUsed.contains(edgeVariables.inverse().get(edgeVariablesAsArray[i])) ? 1 : 0;
}
return edgeVals;
}
private Map<E,Double> getNonZeroEdgeWeights(double[] edgeValues){
Map<E,Double> edgeWeights = new HashMap<E,Double>();
for(int i = 0; i < edgeValues.length; i++){
if(edgeValues[i] > Util.epsilon){
edgeWeights.put(edgeVariables.inverse().get(edgeVariablesAsArray[i]), edgeValues[i]);
}
}
return edgeWeights;
}
private class IntegerCutCallback extends LazyConstraintCallback{
private double minCutVal;
public IntegerCutCallback(){
minCutVal = 1.99;
}
@Override
protected void main() throws IloException {
Set<E> edgesUsed = edgesUsed(this.getValues(edgeVariablesAsArray));
Set<Set<V>> connectedComponents = tspInstance.getConnectedComponents(edgesUsed);
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, tspInstance.cutEdges(connectedComponent)),2));
}
}
}
}
private class MinCutCallback extends UserCutCallback{
private BackOffFunction backOff;
private double minCutVal;
public MinCutCallback(){
minCutVal = 1.99;
this.backOff = new BackOffFunction.QuadraticWaiting();
}
@Override
protected void main() throws IloException {
//do not attempt to add any cutset constraints until CPLEX is done adding constraints
if(!this.isAfterCutLoop()){
return;
}
double[] edgeVals = this.getValues(edgeVariablesAsArray);
Map<E,Double> edgeWeights = getNonZeroEdgeWeights(edgeVals);
UndirectedGraph<V,E> graph = tspInstance.getGraph();
MinCutSolver<V,E> minCutSolver = new MinCutSolver<V,E>(graph,edgeWeights);
Set<Set<E>> cutSets = new HashSet<Set<E>>();
Set<Set<V>> connectedComponents = minCutSolver.checkConnectedComponents();
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
cutSets.add(tspInstance.cutEdges(connectedComponent));
}
System.out.println("trivial cuts exist, found " + cutSets.size());
}
else if(this.getNnodes64() == 0 || this.backOff.makeAttempt()){
V source = graph.getVertices().iterator().next();
double best = Double.MAX_VALUE;
System.out.print("looking for cut... ");
for(V sink : graph.getVertices()){
if(sink != source){
Set<V> cutNodes = minCutSolver.computeMinCut(source, sink);
Set<E> cutEdges = tspInstance.cutEdges(cutNodes);
double cutVal = Util.calcSum(cutEdges,edgeWeights);
best = Math.min(best, cutVal);
if(cutVal < minCutVal){
cutSets.add(cutEdges);
}
}
}
System.out.println("found " + cutSets.size() + " cuts, best " + best);
}
for(Set<E> cutEdges: cutSets){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, cutEdges), 2));
}
}
}
private class RandomMinCutCallback extends UserCutCallback{
private int numAttempts;
private Random rand;
private BackOffFunction backOff;
private double minCutVal;
public RandomMinCutCallback(){
numAttempts = 20;
rand = new Random();
backOff = new BackOffFunction.QuadraticWaiting();
minCutVal = 1.99;
}
@Override
protected void main() throws IloException {
//do not attempt to add any cutset constraints until CPLEX is done adding constraints
if(!this.isAfterCutLoop()){
return;
}
double[] edgeVals = this.getValues(edgeVariablesAsArray);
Map<E,Double> edgeWeights = getNonZeroEdgeWeights(edgeVals);
UndirectedGraph<V,E> graph = tspInstance.getGraph();
MinCutSolver<V,E> minCutSolver = new MinCutSolver<V,E>(graph,edgeWeights);
Set<Set<E>> cutSets = new HashSet<Set<E>>();
Set<Set<V>> connectedComponents = minCutSolver.checkConnectedComponents();
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
cutSets.add(tspInstance.cutEdges(connectedComponent));
}
System.out.println("trivial cuts exist, found " + cutSets.size());
}
else if(this.getNnodes64() == 0 || this.backOff.makeAttempt()){
double best = Double.MAX_VALUE;
System.out.print("looking for cut... ");
List<V> nodes = new ArrayList<V>(graph.getVertices());
for(int i = 0; i <numAttempts; i++){
V source = nodes.get(rand.nextInt(nodes.size()));
V sink = nodes.get(rand.nextInt(nodes.size()));
if(source != sink){
Set<V> cutNodes = minCutSolver.computeMinCut(source, sink);
Set<E> cutEdges = tspInstance.cutEdges(cutNodes);
double cutVal = Util.calcSum(cutEdges,edgeWeights);
best = Math.min(best, cutVal);
if(cutVal < minCutVal){
cutSets.add(cutEdges);
}
}
}
System.out.println("found " + cutSets.size() + " cuts, best " + best);
}
for(Set<E> cutEdges: cutSets){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, cutEdges), 2));
}
}
}
}