Inventory Routing (IRP)¶

using System;
using System.IO;
using System.Globalization;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

    // Horizon length
    int horizonLength;

    // Capacity
    int capacity;

    // Start level at the supplier
    int startLevelSupplier;

    // Production rate of the supplier
    int productionRateSupplier;

    // Holding costs of the supplier
    double holdingCostSupplier;

    // Start level of the customers
    int[] startLevel;

    // Max level of the customers
    int[] maxLevel;

    // Demand rate of the customers
    int[] demandRate;

    // Holding costs of the customers
    double[] holdingCost;

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

    // Distances between customers and supplier
    long[] distSupplierData;

    // Decision variables
    private HxExpression[][] delivery;

    // Decision variables
    private HxExpression[] route;

    // Are the customers receiving products
    private HxExpression[][] isDelivered;

    // Distance traveled by each truck
    HxExpression[] distRoutes;

    // Inventory at the supplier
    HxExpression[] inventorySupplier;

    // Inventory at the customers
    HxExpression[][] inventory;

    // Total inventory cost at the supplier
    HxExpression totalCostInventorySupplier;

    // Inventory cost at a customer
    HxExpression[] costInventoryCustomer;

    // Total inventory cost at customers
    HxExpression totalCostInventory;

    // Total transportation cost
    HxExpression totalCostRoute;

    // Objective
    HxExpression objective;

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

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

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

        delivery = new HxExpression[horizonLength][];
        route = new HxExpression[horizonLength];
        isDelivered = new HxExpression[horizonLength][];
        distRoutes = new HxExpression[horizonLength];

        // Quantity of product delivered at each discrete time instant of
        // the planning time horizon to each customer
        for (int t = 0; t < horizonLength; ++t)
        {
            delivery[t] = new HxExpression[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
                delivery[t][i] = model.Float(0, capacity);
        }

        // Sequence of customers visited at each discrete time instant of
        // the planning time horizon
        for (int t = 0; t < horizonLength; ++t)
            route[t] = model.List(nbCustomers);

        // Create distances as arrays to be able to access them with an "at" operator
        HxExpression distSupplier = model.Array(distSupplierData);
        HxExpression distMatrix = model.Array(distMatrixData);

        for (int t = 0; t < horizonLength; ++t)
        {
            HxExpression sequence = route[t];
            HxExpression c = model.Count(sequence);

            isDelivered[t] = new HxExpression[nbCustomers];
            // Customers receive products only if they are visited
            for (int i = 0; i < nbCustomers; ++i)
                isDelivered[t][i] = model.Contains(sequence, i);

            // Distance traveled at instant t
            HxExpression distLambda = model.LambdaFunction(
                i => distMatrix[sequence[i - 1], sequence[i]]
            );
            distRoutes[t] = model.If(
                c > 0,
                distSupplier[sequence[0]]
                    + model.Sum(model.Range(1, c), distLambda)
                    + distSupplier[sequence[c - 1]],
                0
            );
        }

        // Stockout constraints at the supplier
        inventorySupplier = new HxExpression[horizonLength + 1];
        inventorySupplier[0] = model.CreateConstant(startLevelSupplier);
        for (int t = 1; t < horizonLength + 1; ++t)
        {
            inventorySupplier[t] =
                inventorySupplier[t - 1] - model.Sum(delivery[t - 1]) + productionRateSupplier;
            if (t != horizonLength)
                model.Constraint(inventorySupplier[t] >= model.Sum(delivery[t]));
        }

        // Stockout constraints at the customers
        inventory = new HxExpression[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
        {
            inventory[i] = new HxExpression[horizonLength + 1];
            inventory[i][0] = model.CreateConstant(startLevel[i]);
            for (int t = 1; t < horizonLength + 1; ++t)
            {
                inventory[i][t] = inventory[i][t - 1] + delivery[t - 1][i] - demandRate[i];
                model.Constraint(inventory[i][t] >= 0);
            }
        }

        for (int t = 0; t < horizonLength; ++t)
        {
            // Capacity constraints
            model.Constraint(model.Sum(delivery[t]) <= capacity);

            // Maximum level constraints
            for (int i = 0; i < nbCustomers; ++i)
            {
                model.Constraint(delivery[t][i] <= maxLevel[i] - inventory[i][t]);
                model.Constraint(delivery[t][i] <= maxLevel[i] * isDelivered[t][i]);
            }
        }

        // Total inventory cost at the supplier
        totalCostInventorySupplier = holdingCostSupplier * model.Sum(inventorySupplier);

        // Total inventory cost at customers
        costInventoryCustomer = new HxExpression[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            costInventoryCustomer[i] = holdingCost[i] * model.Sum(inventory[i]);
        totalCostInventory = model.Sum(costInventoryCustomer);

        // Total transportation cost
        totalCostRoute = model.Sum(distRoutes);

        // Objective: minimize the sum of all costs
        objective = model.Sum(
            model.Sum(totalCostInventorySupplier, totalCostInventory),
            totalCostRoute
        );
        model.Minimize(objective);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     * - total distance run by the vehicle
     * - the nodes visited at each time step (omitting the start/end at the supplier) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalCostRoute.GetValue());
            for (int t = 0; t < horizonLength; ++t)
            {
                HxCollection routeCollection = route[t].GetCollectionValue();
                for (int i = 0; i < routeCollection.Count(); ++i)
                    output.Write((routeCollection[i] + 1) + " ");
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Irp 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 (Irp model = new Irp())
        {
            model.ReadInputIrp(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

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

    // The input files follow the "Archetti" format
    private void ReadInputIrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            splitted = ReadNextLine(input);
            nbCustomers = int.Parse(splitted[0]) - 1;
            horizonLength = int.Parse(splitted[1]);
            capacity = int.Parse(splitted[2]);
            splitted = ReadNextLine(input);
            Console.WriteLine(1);
            Console.WriteLine(splitted[1]);
            double xCoordSupplier = double.Parse(splitted[1], CultureInfo.InvariantCulture);
            Console.WriteLine(2);
            double yCoordSupplier = double.Parse(splitted[2], CultureInfo.InvariantCulture);
            Console.WriteLine(3);
            startLevelSupplier = int.Parse(splitted[3]);
            Console.WriteLine(4);
            productionRateSupplier = int.Parse(splitted[4]);
            holdingCostSupplier = double.Parse(splitted[5], CultureInfo.InvariantCulture);
            double[] xCoord = new double[nbCustomers];
            double[] yCoord = new double[nbCustomers];

            startLevel = new int[nbCustomers];
            maxLevel = new int[nbCustomers];
            demandRate = new int[nbCustomers];
            holdingCost = new double[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
            {
                splitted = ReadNextLine(input);
                xCoord[i] = double.Parse(splitted[1], CultureInfo.InvariantCulture);
                yCoord[i] = double.Parse(splitted[2], CultureInfo.InvariantCulture);
                startLevel[i] = int.Parse(splitted[3]);
                maxLevel[i] = int.Parse(splitted[4]);
                demandRate[i] = int.Parse(splitted[6]);
                holdingCost[i] = double.Parse(splitted[7], CultureInfo.InvariantCulture);
            }

            ComputeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);
        }
    }

    // Compute the distance matrix
    private void ComputeDistanceMatrix(
        double xCoordSupplier,
        double yCoordSupplier,
        double[] xCoord,
        double[] yCoord
    )
    {
        distMatrixData = new long[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new long[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                long dist = ComputeDist(xCoord[i], xCoord[j], yCoord[i], yCoord[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distSupplierData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            distSupplierData[i] = ComputeDist(xCoordSupplier, xCoord[i], yCoordSupplier, yCoord[i]);
    }

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