Inventory Routing (IRP)¶

/********** Irp.cs **********/

using System;
using System.IO;
using System.Globalization;
using localsolver;

public class Irp : IDisposable
{
    // Solver
    LocalSolver localsolver;

    // 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[][] distanceMatrix;

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

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

    // Decision variables
    private LSExpression[] route;

    // Are the customers receiving products
    private LSExpression[][] isDelivered;
    
    // Distance traveled by each truck
    LSExpression[] costRoute;

    // Inventory at the supplier
    LSExpression[] inventorySupplier;

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

    // Total inventory cost at the supplier
    LSExpression totalCostInventorySupplier;

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

    // Total inventory cost at customers
    LSExpression totalCostInventory;

    // Total transportation cost
    LSExpression totalCostRoute;

    // Objective
    LSExpression objective;

    public Irp()
    {
        localsolver = new LocalSolver();
    }

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

    void Solve(int limit)
    {
        // Declares the optimization model.
        LSModel model = localsolver.GetModel();

        delivery = new LSExpression[horizonLength][];
        route =  new LSExpression[horizonLength];
        isDelivered = new LSExpression[horizonLength][];
        costRoute = new LSExpression[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 LSExpression[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 it with an "at" operator
        LSExpression distanceSupplierArray = model.Array(distanceSupplier);
        LSExpression distanceArray = model.Array(distanceMatrix);

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

            isDelivered[t] = new LSExpression[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
            LSExpression distSelector = model.LambdaFunction(i => distanceArray[sequence[i - 1], sequence[i]]);
            costRoute[t] = model.If(c > 0, distanceSupplierArray[sequence[0]] + model.Sum(model.Range(1, c), distSelector) + distanceSupplierArray[sequence[c - 1]], 0);
        }

        // Stockout constraints at the supplier
        inventorySupplier =  new LSExpression[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 LSExpression[nbCustomers][];
        for (int i = 0; i < nbCustomers; i++){
            inventory[i] = new LSExpression[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 LSExpression[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(costRoute);
        
        // Objective: minimize the sum of all costs
        objective = model.Sum(model.Sum(totalCostInventorySupplier, totalCostInventory), totalCostRoute);
        model.Minimize(objective);

        model.Close();

        // Parameterizes the solver.
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.Solve();
    }

    // Writes 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++)
            {
                LSCollection 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);
        }
    }

    // The input files follow the "Archetti" format.
    private void ReadInputIrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            splitted = input.ReadLine().Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
            nbCustomers = int.Parse(splitted[0]) - 1;
            horizonLength = int.Parse(splitted[1]);
            capacity = int.Parse(splitted[2]);
            splitted = input.ReadLine().Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
            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 = input.ReadLine().Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
                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);
            }
            
            // Compute distance matrix
            ComputeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);
        }
    }

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

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

        distanceSupplier = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distanceSupplier[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))));
    }
}