Vehicle Routing Problem with Transhipment Facilities (VRPTF)¶

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

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

    // Number of customers
    int nbCustomers;

    // Customers coordinates
    int[] xCustomers;
    int[] yCustomers;

    // Customers demands
    int[] demandsData;

    // Number of facilities
    int nbFacilities;

    // Depots coordinates
    int[] xFacilities;
    int[] yFacilities;

    // Number of points
    int nbPoints;

    // Depot coordinates
    int xDepot;
    int yDepot;

    // Number of trucks
    int nbTrucks;

    // Capacity of a truck
    int truckCapacity;

    // Distance matrix
    long[][] distMatrixData;

    // Distance to depot array
    long[] distDepotsData;

    // Assignement costs
    long[] assignmentCostsData;

    // Decision variables
    HxExpression[] routesSequences;

    // Objective value
    HxExpression totalDistanceCost;
    HxExpression totalAssignementCost;
    HxExpression totalCost;

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

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

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

        double totalDemand = 0;
        for (int c = 0; c < nbCustomers; ++c)
        {
            totalDemand += demandsData[c];
        }
        int minNbTrucks = (int)Math.Ceiling(totalDemand / truckCapacity);
        nbTrucks = (int)Math.Ceiling(1.5 * minNbTrucks);

        routesSequences = new HxExpression[nbTrucks];

        // A route is represented as a list containing the points in the order they are visited
        for (int r = 0; r < nbTrucks; ++r)
        {
            routesSequences[r] = m.List(nbPoints);
        }

        HxExpression routes = m.Array(routesSequences);

        // Each point must be visited at most once
        m.Constraint(m.Disjoint(routesSequences));

        HxExpression demands = m.Array(demandsData);
        HxExpression distMatrix = m.Array(distMatrixData);
        HxExpression distDepots = m.Array(distDepotsData);
        HxExpression assignmentCosts = m.Array(assignmentCostsData);

        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] assignmentCostRoutes = new HxExpression[nbTrucks];

        for (int c = 0; c < nbCustomers; ++c)
        {
            int startFacilities = nbCustomers + c * nbFacilities;
            int endFacilities = startFacilities + nbFacilities;

            // Each customer is either contained in a route or assigned to a facility
            HxExpression facilityUsedSum = m.Sum();
            for (int f = startFacilities; f < endFacilities; ++f)
            {
                facilityUsedSum.AddOperand(m.Contains(routes, f));
            }
            m.Constraint(m.Contains(routes, c) + facilityUsedSum == 1);
        }

        for (int r = 0; r < nbTrucks; ++r)
        {
            HxExpression route = routesSequences[r];
            HxExpression c = m.Count(route);

            // Each truck cannot carry more than its capacity
            HxExpression demandLambda = m.LambdaFunction(j => m.At(demands, j));
            HxExpression quantityServed = m.Sum(route, demandLambda);
            m.Constraint(quantityServed <= truckCapacity);

            HxExpression distLambda = m.LambdaFunction(i => m.At(distMatrix, m.At(route, i), m.At(route, i + 1)));

            // Truck is used if it visits at least one point
            HxExpression truckUsed = c > 0;

            // Distance traveled by each truck
            distRoutes[r] = m.Sum(m.Range(0, c - 1), distLambda) +
                    m.If(truckUsed,
                    m.At(distDepots, m.At(route, 0)) +
                            m.At(distDepots, m.At(route, c - 1)),
                            0);

            // The cost to assign customers to their facility
            HxExpression assignmentCostLambda = m.LambdaFunction(i => m.At(assignmentCosts, i));
            assignmentCostRoutes[r] = m.Sum(route, assignmentCostLambda);
        }

        // The total distance travelled
        totalDistanceCost = m.Sum(distRoutes);
        // The total assignement cost
        totalAssignementCost = m.Sum(assignmentCostRoutes);

        // Objective: minimize the sum of the total distance travelled and the total assignement cost
        totalCost = totalDistanceCost + totalAssignementCost;

        m.Minimize(totalCost);
        m.Close();

        optimizer.GetParam().SetTimeLimit(limit);
        optimizer.Solve();
    }

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] line;
            line = ReadNextLine(input);
            nbCustomers = int.Parse(line[0]);
            xCustomers = new int[nbCustomers];
            yCustomers = new int[nbCustomers];

            line = ReadNextLine(input);
            nbFacilities = int.Parse(line[0]);
            xFacilities = new int[nbFacilities];
            yFacilities = new int[nbFacilities];

            // A point is either a customer or a facility
            // Facilities are duplicated for each customer
            nbPoints = nbCustomers + nbCustomers * nbFacilities;
            distMatrixData = new long[nbPoints][];
            distDepotsData = new long[nbPoints];
            demandsData = new int[nbPoints];

            line = ReadNextLine(input);
            for (int f = 0; f < nbFacilities; ++f)
            {
                line = ReadNextLine(input);
                xFacilities[f] = int.Parse(line[0]);
                yFacilities[f] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            for (int c = 0; c < nbCustomers; ++c)
            {
                line = ReadNextLine(input);
                xCustomers[c] = int.Parse(line[0]);
                yCustomers[c] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            truckCapacity = int.Parse(line[0]);
            line = ReadNextLine(input);

            // Facility capacities : skip
            for (int f = 0; f < nbFacilities; ++f)
            {
                ReadNextLine(input);
            }
            line = ReadNextLine(input);
            for (int c = 0; c < nbCustomers; ++c)
            {
                line = ReadNextLine(input);
                demandsData[c] = int.Parse(line[0]);
                for (int f = 0; f < nbFacilities; f++)
                {
                    demandsData[nbCustomers + c * nbFacilities + f] = demandsData[c];
                }
            }
            line = ReadNextLine(input);

            ComputeDepotCoordinates();
            ComputeDistanceArrays();
            ComputeAssignmentCosts();
        }
    }

    void ComputeDepotCoordinates()
    {
        // Compute the coordinates of the bounding box containing all of the points
        int xMin = Math.Min(xCustomers.Min(), xFacilities.Min());
        int xMax = Math.Max(xCustomers.Max(), xFacilities.Max());
        int yMin = Math.Min(yCustomers.Min(), yFacilities.Min());
        int yMax = Math.Max(yCustomers.Max(), yFacilities.Max());

        // We assume that the depot is at the center of the bounding box
        xDepot = xMin + (xMax - xMin) / 2;
        yDepot = yMin + (yMax - yMin) / 2;
    }

    void ComputeDistanceArrays()
    {
        distDepotsData = new long[nbPoints];

        // Customer to depot
        for (int c = 0; c < nbCustomers; ++c)
        {
            distDepotsData[c] = ComputeDist(xCustomers[c], xDepot, yCustomers[c], yDepot);
        }

        // Facility to depot
        for (int c = 0; c < nbCustomers; ++c)
        {
            for (int f = 0; f < nbFacilities; ++f)
            {
                distDepotsData[nbCustomers + c * nbFacilities + f] = ComputeDist(xFacilities[f], xDepot, yFacilities[f], yDepot);
            }
        }

        distMatrixData = new long[nbPoints][];
        for (int i = 0; i < nbPoints; i++)
        {
            distMatrixData[i] = new long[nbPoints];
        }

        // Distances between customers
        for (int c1 = 0; c1 < nbCustomers; ++c1)
        {
            for (int c2 = 0; c2 < nbCustomers; ++c2)
            {
                long dist = ComputeDist(xCustomers[c1], xCustomers[c2], yCustomers[c1], yCustomers[c2]);
                distMatrixData[c1][c2] = dist;
                distMatrixData[c2][c1] = dist;
            }
        }

        // Distances between customers and facilities
        for (int c1 = 0; c1 < nbCustomers; ++c1)
        {
            for (int f = 0; f < nbFacilities; ++f)
            {
                long dist = ComputeDist(xFacilities[f], xCustomers[c1], yFacilities[f], yCustomers[c1]);
                for (int c2 = 0; c2 < nbCustomers; ++c2)
                {
                    // Index representing serving c2 through facility f
                    int facilityIndex = nbCustomers + c2 * nbFacilities + f;
                    distMatrixData[facilityIndex][c1] = dist;
                    distMatrixData[c1][facilityIndex] = dist;
                }
            }
        }
    
        // Distances between facilities
        for (int f1 = 0; f1 < nbFacilities; ++f1)
        {
            for (int f2 = 0; f2 < nbFacilities; ++f2)
            {
                long dist = ComputeDist(xFacilities[f1], xFacilities[f2], yFacilities[f1], yFacilities[f2]);
                for (int c1 = 0; c1 < nbCustomers; ++c1)
                {
                    // Index representing serving c1 through facility f1
                    int index1 = nbCustomers + c1 * nbFacilities + f1;
                    for (int c2 = 0; c2 < nbCustomers; ++c2)
                    {
                        // Index representing serving c2 through facility f2
                        int index2 = nbCustomers + c2 * nbFacilities + f2;
                        distMatrixData[index1][index2] = dist;
                    }
                }
            }
        }
    }

    private long ComputeDist(int xi, int xj, int yi, int yj)
    {
        double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        return Convert.ToInt64(Math.Round(exactDist));
    }

    private void ComputeAssignmentCosts() {
        // Compute assignment cost for each point
        assignmentCostsData = new long[nbPoints];
        for (int c = 0; c < nbCustomers; ++c)
        {
            assignmentCostsData[c] = 0;
            for (int f = 0; f < nbFacilities; ++f)
            {
                // Cost of serving customer c through facility f
                assignmentCostsData[nbCustomers + c * nbFacilities + f] =
                    distMatrixData[c][nbCustomers + c * nbFacilities + f];
            }
        }
    }

    /* Write the solution in a file */
    void WriteSolution(string infile, string outfile)
    {
        using (StreamWriter output = new StreamWriter(outfile))
        {
            output.WriteLine(
                "File name: " + infile + "; totalCost = " + totalCost.GetIntValue()
                + "; totalDistance = " + totalDistanceCost.GetIntValue()
                + "; totalAssignementCost = " + totalAssignementCost.GetIntValue()
            );
            for (int r = 0; r < nbTrucks; ++r)
            {
                HxCollection route = routesSequences[r].GetCollectionValue();
                if (route.Count() == 0) continue;
                output.Write("Route " + r + " [");
                for (int i = 0; i < route.Count(); i++)
                {
                    long point = route.Get(i);
                    if (point < nbCustomers)
                    {
                        output.Write("Customer " + point);
                    }
                    else
                    {
                        output.Write("Facility " + point % nbCustomers + " assigned to Customer " + (point - nbCustomers) / nbFacilities);
                    }
                    if (i < route.Count() - 1)
                    {
                        output.Write(", ");
                    }
                }
                output.WriteLine("]");
            }
        }
    }

    String[] ReadNextLine(StreamReader input)
    {
        return input.ReadLine().Split(new[] { '\t' }, StringSplitOptions.RemoveEmptyEntries);
    }

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

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