Vehicle Routing with Time Windows (CVRPTW)¶

using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

public class Cvrptw : IDisposable
{
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Latest allowed arrival to depot
    int maxHorizon;

    // Demand for each customer
    List<int> demandsData;

    // Earliest arrival for each customer
    List<int> earliestStartData;

    // Latest departure from each customer
    List<int> latestEndData;

    // Service time for each customer
    List<int> serviceTimeData;

    // Distance matrix between customers
    double[][] distMatrixData;

    // Distances between customers and depot
    double[] distDepotData;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    HxExpression[] customersSequences;

    // Are the trucks actually used
    HxExpression[] trucksUsed;

    // Cumulated lateness in the solution (must be 0 for the solution to be valid)
    HxExpression totalLateness;

    // Number of trucks used in the solution
    HxExpression nbTrucksUsed;

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    public Cvrptw()
    {
        optimizer = new HexalyOptimizer();
    }

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputCvrptw(fileName);
    }

    public void Dispose()
    {
        if (optimizer != null)
            optimizer.Dispose();
    }

    void Solve(int limit)
    {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        trucksUsed = new HxExpression[nbTrucks];
        customersSequences = new HxExpression[nbTrucks];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] endTime = new HxExpression[nbTrucks];
        HxExpression[] homeLateness = new HxExpression[nbTrucks];
        HxExpression[] lateness = new HxExpression[nbTrucks];

        // Sequence of customers visited by each truck
        for (int k = 0; k < nbTrucks; ++k)
            customersSequences[k] = model.List(nbCustomers);

        // All customers must be visited by exactly one truck
        model.Constraint(model.Partition(customersSequences));

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.Array(demandsData);
        HxExpression earliest = model.Array(earliestStartData);
        HxExpression latest = model.Array(latestEndData);
        HxExpression serviceTime = model.Array(serviceTimeData);
        HxExpression distDepot = model.Array(distDepotData);
        HxExpression distMatrix = model.Array(distMatrixData);

        for (int k = 0; k < nbTrucks; ++k)
        {
            HxExpression sequence = customersSequences[k];
            HxExpression c = model.Count(sequence);

            // A truck is used if it visits at least one customer
            trucksUsed[k] = c > 0;

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression demandLambda = model.LambdaFunction(j => demands[j]);
            HxExpression routeQuantity = model.Sum(sequence, demandLambda);
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            HxExpression distLambda = model.LambdaFunction(
                i => distMatrix[sequence[i - 1], sequence[i]]
            );
            distRoutes[k] =
                model.Sum(model.Range(1, c), distLambda)
                + model.If(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

            // End of each visit
            HxExpression endTimeLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        earliest[sequence[i]],
                        model.If(
                            i == 0,
                            distDepot[sequence[0]],
                            prev + distMatrix[sequence[i - 1], sequence[i]]
                        )
                    ) + serviceTime[sequence[i]]
            );

            endTime[k] = model.Array(model.Range(0, c), endTimeLambda, 0);

            // Arriving home after max_horizon
            homeLateness[k] = model.If(
                trucksUsed[k],
                model.Max(0, endTime[k][c - 1] + distDepot[sequence[c - 1]] - maxHorizon),
                0
            );

            // Completing visit after latest_end
            HxExpression lateLambda = model.LambdaFunction(
                i => model.Max(endTime[k][i] - latest[sequence[i]], 0)
            );
            lateness[k] = homeLateness[k] + model.Sum(model.Range(0, c), lateLambda);
        }

        // Total lateness
        totalLateness = model.Sum(lateness);

        // Total number of trucks used
        nbTrucksUsed = model.Sum(trucksUsed);

        // Total distance traveled (convention in Solomon's instances is to round to 2 decimals)
        totalDistance = model.Round(100 * model.Sum(distRoutes)) / 100;

        // Objective: minimize the lateness, then the number of trucks used, then the distance traveled
        model.Minimize(totalLateness);
        model.Minimize(nbTrucksUsed);
        model.Minimize(totalDistance);

        model.Close();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - number of trucks used and total distance
     *  - for each truck the customers visited (omitting the start/end at the depot) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(nbTrucksUsed.GetValue() + " " + totalDistance.GetDoubleValue());
            for (int k = 0; k < nbTrucks; ++k)
            {
                if (trucksUsed[k].GetValue() != 1)
                    continue;
                // Values in sequence are in 0...nbCustomers. +1 is to put it back in 1...nbCustomers+1
                // as in the data files (0 being the depot)
                HxCollection customersCollection = customersSequences[k].GetCollectionValue();
                for (int i = 0; i < customersCollection.Count(); ++i)
                    output.Write((customersCollection[i] + 1) + " ");
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Cvrptw inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

        using (Cvrptw model = new Cvrptw())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

    private string[] SplitInput(StreamReader input)
    {
        string line = input.ReadLine();
        if (line == null)
            return new string[0];
        return line.Split(new[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
    }

    // The input files follow the "Solomon" format
    private void ReadInputCvrptw(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;

            input.ReadLine();
            input.ReadLine();
            input.ReadLine();
            input.ReadLine();

            splitted = SplitInput(input);
            nbTrucks = int.Parse(splitted[0]);
            truckCapacity = int.Parse(splitted[1]);

            input.ReadLine();
            input.ReadLine();
            input.ReadLine();
            input.ReadLine();

            splitted = SplitInput(input);
            int depotX = int.Parse(splitted[1]);
            int depotY = int.Parse(splitted[2]);
            maxHorizon = int.Parse(splitted[5]);

            List<int> customersX = new List<int>();
            List<int> customersY = new List<int>();
            demandsData = new List<int>();
            earliestStartData = new List<int>();
            latestEndData = new List<int>();
            serviceTimeData = new List<int>();

            while (!input.EndOfStream)
            {
                splitted = SplitInput(input);
                if (splitted.Length < 7)
                    break;
                customersX.Add(int.Parse(splitted[1]));
                customersY.Add(int.Parse(splitted[2]));
                demandsData.Add(int.Parse(splitted[3]));
                int ready = int.Parse(splitted[4]);
                int due = int.Parse(splitted[5]);
                int service = int.Parse(splitted[6]);

                earliestStartData.Add(ready);
                latestEndData.Add(due + service); // in input files due date is meant as latest start time
                serviceTimeData.Add(service);
            }

            nbCustomers = customersX.Count;

            ComputeDistanceMatrix(depotX, depotY, customersX, customersY);
        }
    }

    // Compute the distance matrix
    private void ComputeDistanceMatrix(
        int depotX,
        int depotY,
        List<int> customersX,
        List<int> customersY
    )
    {
        distMatrixData = new double[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new double[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                double dist = ComputeDist(
                    customersX[i],
                    customersX[j],
                    customersY[i],
                    customersY[j]
                );
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distDepotData = new double[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            distDepotData[i] = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
    }

    private double ComputeDist(int xi, int xj, int yi, int yj)
    {
        return Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
    }
}