Exact Optimization

Google OR-Tools

According to Google, "​OR-Tools is an open source software suite for optimization, tuned for tackling the world's toughest problems in vehicle routing, flows, integer and linear programming, and constraint programming." OR-Tools are available for C++, Python, Java and .NET, the latter API can be used in HeuristicLab.

OR-Tools includes solvers for

• Contraint Programming (​CP, ​CP-SAT),
• Linear and Mixed-Integer Programming (​BOP, ​Cbc, ​Clp, ​CPLEX, ​Glop, ​GLPK, ​Gurobi and ​SCIP),
• Routing,
• Scheduling and
• Network Flows.

For a complete overview, see ​https://developers.google.com/optimization/introduction/overview. You can use all solvers in HeuristicLab, just create a C# Script and make sure the ​Microsoft Visual C++ 2017 Redistributable (x64) is installed. Note that most linear and mixed-integer solvers must be installed (and licensed) separately, see the next section for details. The HeuristicLab.OrTools plugin is currently only available on 64-bit Windows.

Mixed-Integer (Linear) Programming (MIP, MILP)

OR-Tools contains a ​linear solver wrapper that can be used to define a MIP model once and solve it with a range of MIP solvers. Some solvers are already built into OR-Tools (BOP, Cbc, Clp, Glop) but for most solvers OR-Tools have to be built from source code, this version then has a dependency on the external solver. HeurisicLab includes a custom version of OR-Tools that allows users to simply install additional solvers without depending on them. Download the solvers here, they are free for academic use:

• ​CPLEX,
• ​Gurobi,
• ​GLPK (currently not supported by HeuristicLab) and
• ​SCIP.

Make sure that the solver you install is added to the environment variable Path and the LibraryName property of the solver in he Algorithm view is set to the correct name of the solver DLL (e.g. cplex1280.dll for version 12.8.0 of CPLEX). The default values for the solver DLLs can be set in the configuration file of HeuristicLab. Note that most solvers offer to change the ProblemType property to LinearProgramming, the solver then solves the linear programming relaxation of your model (all integer constraints are removed).

For a general introduction to MIP and a simple model implemented using the linear solver wrapper, see ​https://developers.google.com/optimization/mip/integer_opt. For further examples demonstrating more advanced usage of the linear solver wrapper using the C# API, see ​https://github.com/google/or-tools/tree/master/examples/dotnet.

The HeuristicLab.ExactOptimization plugin provides a GUI for MIP, just create a new Mixed-Integer Linear Programming (LP, MIP) algorithm, found under Algorithms, Exact. You can either load your model from a file (MPS, Google OR-Tools Protocol Buffers files) or define your model using C# and directly accessing the linear solver wrapper. The next section shows the implementation of a Knapsack Problem (KSP) and a Traveling Salesman Problem (TSP).

Defining a MIP Model

A programmable MIP model requires you to implement the BuildModel method where you must define the decision variables, the constraints and the objective of your model. The Analyze method is used to retrieve the solution values of the decision variables, which can be added to the results to be accessible in the GUI. The Initialize method can be used for initialization.

Knapsack Problem (KSP)

The following MIP model can be pasted in HeuristicLab in the Problem tab of a Mixed-Integer Linear Programming (LP, MIP) algorithm. The Analyze method shows the definition of the ​0-1 Knapsack Problem. Note that you can drag any instance of a Knapsack Problem (KSP) to the variables list and rename it to Problem (remove the default problem first).

using System;
using System.Linq;
using System.Text;
using Google.OrTools.LinearSolver;
using HeuristicLab.Data;
using HeuristicLab.Encodings.BinaryVectorEncoding;
using HeuristicLab.ExactOptimization.LinearProgramming;
using HeuristicLab.Optimization;
using HeuristicLab.Problems.Knapsack;

public class KnapsackLinearModel : CompiledProblemDefinition, ILinearProblemDefinition {
  private Variable[] x;

  public override void Initialize() {
    if (vars.Contains("Problem"))
      return;
    vars["Problem"] = new KnapsackProblem();
  }

  public void BuildModel(Solver solver) {
    // Retrieve the problem data
    var W = vars.Problem.KnapsackCapacity.Value;
    IntArray weights = vars.Problem.Weights;
    IntArray values = vars.Problem.Values;
    // Define the decision variables
    x = solver.MakeBoolVarArray(values.Count());
    // Define the constraints
    solver.Add(weights.Select((w, i) => w * x[i]).ToArray().Sum() <= W);
    // Define the objective
    solver.Maximize(values.Select((v, i) => v * x[i]).ToArray().Sum());
  }

  public void Analyze(Solver solver, ResultCollection results) {
    // Retrieve the problem data
    var capacity = vars.Problem.KnapsackCapacity;
    var weights = vars.Problem.Weights;
    var values = vars.Problem.Values;
    // Retrieve the solution values of the objective and the decision variables
    var solution = new BinaryVector(x.Select(xi => xi.SolutionValue() == 1).ToArray());
    var quality = new DoubleValue(solver.Objective().Value());
    // Update the problem
    if (vars.Problem.BestKnownQuality == null || quality.Value > vars.Problem.BestKnownQuality.Value) { 
      vars.Problem.BestKnownSolution = solution;
      vars.Problem.BestKnownQuality = quality;
    }
    // Update the result
    results.AddOrUpdateResult("BestKnapsackSolution", new KnapsackSolution(solution, quality, capacity, weights, values));
  }
}

Traveling Salesman Problem (TSP)

The following MIP model can be pasted in HeuristicLab in the Problem tab of a Mixed-Integer Linear Programming (LP, MIP) algorithm. The Analyze method shows the definition of the ​Miller-Tucker-Zemlin formulation of the Traveling Salesman Problem. Note that you can drag any instance of a Traveling Salesman Problem (TSP) to the variables list and rename it to Problem (remove the default problem first).

using System;
using System.Collections.Generic;
using System.Linq;
using Google.OrTools.LinearSolver;
using HeuristicLab.Data;
using HeuristicLab.Encodings.PermutationEncoding;
using HeuristicLab.ExactOptimization.LinearProgramming;
using HeuristicLab.Optimization;
using HeuristicLab.Problems.TravelingSalesman;

public class TravelingSalesmanLinearModel : CompiledProblemDefinition, ILinearProblemDefinition {
  private Variable[,] x;

  public override void Initialize() {
    if (vars.Contains("Problem"))
      return;
    vars["Problem"] = new TravelingSalesmanProblem();
  }

  public void BuildModel(Solver solver) {
    // Miller-Tucker-Zemlin formulation
    TravelingSalesmanProblem problem = vars.Problem;
    DistanceMatrixHelper.CalculateDistanceMatrix(problem);
    var c = problem.DistanceMatrix;
    var n = c.Rows;
    var N = Enumerable.Range(0, n).ToList();
    // Define the decision variables
    x = solver.MakeBoolVarArray(n, n);
    var u = solver.MakeNumVarArray(n, 0, n - 1);
    // Define the constraints
    foreach (var j in N) {
      solver.Add(N.Where(i => i != j).Select(i => x[i, j]).ToArray().Sum() == 1);
    }
    foreach (var i in N) {
      solver.Add(N.Where(j => j != i).Select(j => x[i, j]).ToArray().Sum() == 1);
    }
    for (var i = 1; i < n; i++) {
      for (var j = 1; j < n; j++) {
        solver.Add(u[i] - u[j] + n * x[i, j] <= n - 1);
      }
    }
    // Define the objective
    solver.Minimize(
      (from i in N
       from j in N
       where i != j
       select c[i, j] * x[i, j]
      ).ToArray().Sum()
    );
  }

  public void Analyze(Solver solver, ResultCollection results) {
    // Retrieve the problem data
    var coordinates = vars.Problem.Coordinates;
    // Retrieve the solution values of the objective and the decision variables
    var solution = new Permutation(PermutationTypes.RelativeUndirected, GetTour(x));
    var quality = new DoubleValue(solver.Objective().Value());
    // Update the problem
    if (vars.Problem.BestKnownQuality == null || quality.Value < vars.Problem.BestKnownQuality.Value) { 
      vars.Problem.BestKnownSolution = solution;
      vars.Problem.BestKnownQuality = quality;
    }
    // Update the result
    results.AddOrUpdateResult("Best TSP Solution", new PathTSPTour(coordinates, solution, quality));
  }
  
  private static int[] GetTour(Variable[,] x) {
    var n = x.GetLength(0);
    var tour = new List<int> { 0 };

    while (true) {
      for (var i = 0; i < n; i++) {
        if (x[tour.Last(), i].SolutionValue() > 0) {
          if (i == tour[0])
            return tour.ToArray();
          tour.Add(i);
          break;
        }
      }
    }
  }

  public static class DistanceMatrixHelper {
    private static double CalculateDistanceEuclideanPath(double x1, double y1, double x2, double y2) {
      return Math.Sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
    }

    private static double CalculateDistanceRoundedEuclideanPath(double x1, double y1, double x2, double y2) {
      return Math.Round(Math.Sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)));
    }

    private static double CalculateDistanceUpperEuclideanPath(double x1, double y1, double x2, double y2) {
      return Math.Ceiling(Math.Sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)));
    }

    private const double PI = 3.141592;
    private const double RADIUS = 6378.388;

    private static double CalculateDistanceGeoPath(double x1, double y1, double x2, double y2) {
      double latitude1, longitude1, latitude2, longitude2;
      double q1, q2, q3;
      double length;

      latitude1 = ConvertToRadian(x1);
      longitude1 = ConvertToRadian(y1);
      latitude2 = ConvertToRadian(x2);
      longitude2 = ConvertToRadian(y2);

      q1 = Math.Cos(longitude1 - longitude2);
      q2 = Math.Cos(latitude1 - latitude2);
      q3 = Math.Cos(latitude1 + latitude2);

      length = (int)(RADIUS * Math.Acos(0.5 * ((1.0 + q1) * q2 - (1.0 - q1) * q3)) + 1.0);
      return (length);
    }

    private static double ConvertToRadian(double x) {
      return PI * (Math.Truncate(x) + 5.0 * (x - Math.Truncate(x)) / 3.0) / 180.0;
    }

    private static double CalculateDistance(ITSPEvaluator evaluator, double x1, double y1, double x2, double y2) {
      if (evaluator is TSPEuclideanPathEvaluator)
        return CalculateDistanceEuclideanPath(x1, y1, x2, y2);
      if (evaluator is TSPRoundedEuclideanPathEvaluator)
        return CalculateDistanceRoundedEuclideanPath(x1, y1, x2, y2);
      if (evaluator is TSPUpperEuclideanPathEvaluator)
        return CalculateDistanceUpperEuclideanPath(x1, y1, x2, y2);
      if (evaluator is TSPGeoPathEvaluator)
        return CalculateDistanceGeoPath(x1, y1, x2, y2);
      throw new InvalidOperationException("Unkown distance measure.");
    }

    public static void CalculateDistanceMatrix(TravelingSalesmanProblem problem) {
      var dm = problem.DistanceMatrix;
      if (dm != null)
        return;
      
      // calculate distance matrix
      var c = problem.Coordinates;
      if (c == null) throw new InvalidOperationException("Neither a distance matrix nor coordinates were given.");
      dm = new DistanceMatrix(c.Rows, c.Rows);
      for (var i = 0; i < dm.Rows; i++) {
        for (var j = 0; j < dm.Columns; j++)      
          dm[i, j] = CalculateDistance(problem.Evaluator, c[i, 0], c[i, 1], c[j, 0], c[j, 1]);
      }
      problem.DistanceMatrix = (DistanceMatrix)dm.AsReadOnly();
    }
  }
}

Implementing New Models in HeuristicLab

The class diagram below shows the simplified architecture of the HeuristicLab.ExactOptimization plugin. To implement a new model that can be easily accessed within HeuristicLab, implement the ILinearProblemDefintion interface. The BuildModel method is used to define your model using the linear solver wrapper of OR-Tools, the Analyze method is used to retrieve the solution values of your decision variables. As this works almost the same as defining a programmable model withing the GUI of HeuristicLab, you can have a look at the examples provided above.