Time Dependent Capacitated Vehicle Routing Problem with Time Windows (TDCVRPTW)¶

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

public class Tdcvrptw : 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;

    // Travel time coefficients for each profile
    static double[] shortDistanceTravelTimeProfile = { 1.00, 2.50, 1.75, 2.50, 1.00 };
    static double[] mediumDistanceTravelTimeProfile = { 1.00, 2.00, 1.50, 2.00, 1.00 };
    static double[] longDistanceTravelTimeProfile = { 1.00, 1.60, 1.10, 1.60, 1.00 };
    static double[][] travelTimeProfileMatrix =
    {
        shortDistanceTravelTimeProfile,
        mediumDistanceTravelTimeProfile,
        longDistanceTravelTimeProfile
    };

    // Distance levels
    static int[] distanceLevels = { 10, 25 };

    // Intervals of the temporal discretization
    static double[] timeIntervalSteps = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };

    // Number of time intervals
    static int nbTimeIntervals = timeIntervalSteps.Length - 1;

    // Number of distances levels
    static int nbDistanceLevels = distanceLevels.Length;

    // Travel time between customers for each day part
    double[][][] travelTimeData;

    // Time interval index for each time unit
    int[] timeToMatrixIdxData;

    // Travel time between customers and depot for each day part
    double[][] travelTimeWarehouseData;

    // 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 Tdcvrptw()
    {
        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);
        HxExpression travelTime = model.Array(travelTimeData);
        HxExpression timeToMatrixIdx = model.Array(timeToMatrixIdxData);
        HxExpression travelTimeWarehouse = model.Array(travelTimeWarehouseData);

        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 according to the traffic
            HxExpression endTimeLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        earliest[sequence[i]],
                        model.If(
                            i == 0,
                            travelTimeWarehouse[sequence[0], timeToMatrixIdx[0]],
                            prev
                                + model.At(
                                    travelTime,
                                    sequence[i - 1],
                                    sequence[i],
                                    timeToMatrixIdx[model.Round(prev)]
                                )
                        )
                    ) + serviceTime[sequence[i]]
            );

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

            // Arriving home after max_horizon
            homeLateness[k] = model.If(
                trucksUsed[k],
                model.Max(
                    0,
                    endTime[k][c - 1]
                        + travelTimeWarehouse[
                            sequence[c - 1],
                            timeToMatrixIdx[model.Round(endTime[k][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();

        // Parameterize 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: Tdcvrptw 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 (Tdcvrptw model = new Tdcvrptw())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

    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
    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
    void ComputeDistanceMatrix(int depotX, int depotY, List<int> customersX, List<int> customersY)
    {
        distMatrixData = new double[nbCustomers][];
        travelTimeData = new double[nbCustomers][][];
        timeToMatrixIdxData = new int[maxHorizon];
        travelTimeWarehouseData = new double[nbCustomers][];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i] = new double[nbCustomers];
            double[][] timeMatrix = new double[nbCustomers][];
            for (int j = 0; j < nbCustomers; ++j)
                timeMatrix[j] = new double[nbTimeIntervals];
            travelTimeData[i] = timeMatrix;
        }

        for (int i = 0; i < nbCustomers; ++i)
            travelTimeWarehouseData[i] = new double[nbTimeIntervals];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int k = 0; k < nbTimeIntervals; ++k)
                travelTimeData[i][i][k] = 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;

                int profileIdx = GetProfile(dist);
                for (int k = 0; k < nbTimeIntervals; ++k)
                {
                    double localTravelTime = travelTimeProfileMatrix[profileIdx][k] * dist;
                    travelTimeData[i][j][k] = localTravelTime;
                    travelTimeData[j][i][k] = localTravelTime;
                }
            }
        }

        for (int i = 0; i < nbTimeIntervals; ++i)
        {
            int timeStepStart = (int)Math.Round(timeIntervalSteps[i] * maxHorizon);
            int timeStepEnd = (int)Math.Round(timeIntervalSteps[i + 1] * maxHorizon);
            for (int j = timeStepStart; j < timeStepEnd; ++j)
                timeToMatrixIdxData[j] = i;
        }

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

            int profileIdx = GetProfile(dist);
            for (int j = 0; j < nbTimeIntervals; ++j)
            {
                double localTravelTimeWarehouse = travelTimeProfileMatrix[profileIdx][j] * dist;
                travelTimeWarehouseData[i][j] = localTravelTimeWarehouse;
            }
        }
    }

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

    int GetProfile(double dist)
    {
        int idx = 0;
        while (idx < nbDistanceLevels && dist > distanceLevels[idx])
            idx += 1;
        return idx;
    }
}