Inventory Routing (IRP)

Principles learned

  • Use a lambda expression to compute a sum with a variable number of terms

  • Use ternary conditions

  • Define a sequence of expressions

  • Access a multi-dimensional array with an “at” operator

Problem

../_images/irp.svg

The inventory routing problem (IRP) is a distribution problem in which a product has to be shipped from a supplier to several customers over a given time horizon. Each customer defines a maximum inventory level. The supplier monitors the inventory of each customer and determines its replenishment policy, guaranteeing that no stockout occurs at the customer (vendor-managed inventory policy). Shipments from the supplier to the customers are performed by a vehicle of given capacity.

Download the example


Data

The instances provided come from the Archetti et al. instances.

The format of the data files is as follows:

  • The first line gives the number of customers, the number of discrete time instants of the planning time horizon and the transportation capacity of the vehicle.

  • The second line contains information concerning the supplier:

    • The index of the supplier

    • The x coordinate

    • The y coordinate

    • The starting level of the inventory

    • The quantity of product made available at each discrete time instant of the planning time horizon

    • The unit inventory cost

  • The next lines contain, for each customer, the following information:

    • The index of the customer

    • The x coordinate

    • The y coordinate

    • The starting level of the inventory

    • The maximum level of the inventory

    • The minimum level of the inventory

    • The quantity absorbed at each discrete time instant of the planning time horizon

    • The unit inventory cost

Program

This LocalSolver model defines a route for each discrete time instant of the planning time horizon t as a list variable (route[t]). It corresponds to the sequence of customers visited.

The inventory level of the supplier at time t is given by the level at time t−1, plus the product quantity made available at time t−1, minus the total quantity shipped to the customers at time t-1. The stockout constraints at the supplier guarantee that for each delivery time t, the inventory level at the supplier is sufficient to ship the total quantity delivered to the customers at time t.

For each customer, the inventory level at time t is given by the level at time t-1, plus the product quantity shipped from the supplier to the customer at time t-1, minus the product quantity consumed at time t-1. The stockout constraints at the customers guarantee that for each customer, the inventory level at each time t is positive.

The capacity constraints guarantee that the total quantity of the product loaded on the vehicle at time t does not exceed the transportation capacity.

The maximum level constraints guarantee that the product quantity shipped from the supplier to each customer is not greater than the maximum inventory level of each customer.

The objective is to minimize the sum of all costs: the total inventory cost at the supplier, the total inventory cost at customers and the total transportation cost.

Execution:
localsolver irp.lsp inFileName=instances/abs3n10.dat [lsTimeLimit=] [solFileName=]
/********** irp.lsp **********/
use io;

/* Reads instance data. The input files follow the "Archetti" format. */
function input() {
    usage = "\nUsage: localsolver irp.lsp " + 
            "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]\n";

    if (inFileName == nil) throw usage;
    
    readInputIrp();

    // Compute distance matrix
    computeDistanceMatrix();
    
}

/* Declares the optimization model. */
function model() {
    // Quantity of product delivered at each discrete time instant of the planning time horizon to each customer
    delivery[t in 0..horizonLength-1][i in 0..nbCustomers-1] <- float(0, capacity);

    // Sequence of customers visited at each discrete time instant of the planning time horizon
    route[0..horizonLength-1] <- list(nbCustomers);

    for [t in 0..horizonLength-1] {
        local sequence <- route[t];
        local c <- count(sequence);

        // Customers receive products only if they are visited
        isDelivered[t][i in 0..nbCustomers-1] <- contains(sequence, i);

        // Distance traveled at instant t
        costRoute[t] <- (c > 0) ? distanceSupplier[sequence[0]] + 
                sum(1..c-1, i => distance[sequence[i-1]][sequence[i]]) + 
                distanceSupplier[sequence[c-1]] : 0;
    }
    
    // Stockout constraints at the supplier
    inventorySupplier[0] <- startLevelSupplier;
    for [t in 1..horizonLength] {
        inventorySupplier[t] <- inventorySupplier[t-1] - sum[i in 0..nbCustomers-1](delivery[t-1][i]) + productionRate_supplier;
        if (t != horizonLength)
            constraint inventorySupplier[t] >= sum[i in 0..nbCustomers-1](delivery[t][i]);
    }

    // Stockout constraints at the customers
    inventory[i in 0..nbCustomers-1][0] <- startLevel[i];
    for [i in 0..nbCustomers-1] {
        for [t in 1..horizonLength] {
            inventory[i][t] <- inventory[i][t-1] + delivery[t-1][i] - demandRate[i];
            constraint inventory[i][t] >= 0;
        }
    }
    
    for [t in 0..horizonLength-1] {
        // Capacity constraints
        constraint sum[i in 0..nbCustomers-1](delivery[t][i]) <= capacity;
        
        // Maximum level constraints
        for [i in 0..nbCustomers-1] {
            constraint delivery[t][i] <= maxLevel[i] - inventory[i][t];
            constraint delivery[t][i] <= maxLevel[i] * isDelivered[t][i];
        }
    }

    // Total inventory cost at the supplier
    costInventorySupplier <- holdingCostSupplier * sum[t in 0..horizonLength](inventorySupplier[t]);

    // Total inventory cost at customers
    costInventory <- sum[i in 0..nbCustomers-1][t in 0..horizonLength](holdingCost[i] * inventory[i][t]);

    // Total transportation cost
    costRoute <- sum[t in 0..horizonLength-1](costRoute[t]);

    // Objective: minimize the sum of all costs
    objective <- costInventorySupplier + costInventory + costRoute;
    minimize(objective);
}

/* Parameterizes the solver. */
function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 20;
}

/* 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) */
function output() {
    if (solFileName == nil) return;
    local outfile = io.openWrite(solFileName);

    outfile.println(costRoute.value);
    for [t in 0..horizonLength-1]{
        for [customer in route[t].value] {
            outfile.print(customer + 1, " ");
        }
        outfile.println();
    }
}

function readInputIrp() {
    inFile = io.openRead(inFileName);

    // General data
    firstLine = inFile.readln().trim().split();
    nbCustomers = firstLine[0].toInt()-1;
    horizonLength = firstLine[1].toInt();
    capacity = firstLine[2].toInt();

    // Supplier data
    secondLine = inFile.readln().trim().split();
    xCoordSupplier = secondLine[1].toDouble();
    yCoordSupplier = secondLine[2].toDouble();
    startLevelSupplier = secondLine[3].toInt();
    productionRate_supplier = secondLine[4].toInt();
    holdingCostSupplier = secondLine[5].toDouble();

    // Customers data
    for[i in 0..nbCustomers-1] {
        elements = inFile.readln().trim().split();
        xCoord[i] = elements[1].toDouble();
        yCoord[i] = elements[2].toDouble();
        startLevel[i] = elements[3].toInt();
        maxLevel[i] = elements[4].toInt();
        minLevel[i] = elements[5].toInt();
        demandRate[i] = elements[6].toInt();
        holdingCost[i] = elements[7].toDouble();
    }
}

function computeDistanceMatrix() {
    distance[i in 0..nbCustomers-1][j in 0..nbCustomers-1] = round(sqrt((xCoord[i] - xCoord[j]) * (xCoord[i] - xCoord[j])
            + (yCoord[i] - yCoord[j]) * (yCoord[i] - yCoord[j])));
    distanceSupplier[i in 0..nbCustomers-1] = round(sqrt((xCoord[i] - xCoordSupplier) * (xCoord[i] - xCoordSupplier)
            + (yCoord[i] - yCoordSupplier) * (yCoord[i] - yCoordSupplier)));
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python irp.py instances\abs3n10.dat
Execution (Linux)
export PYTHONPATH=/opt/localsolver_11_0/bin/python
python irp.py instances/abs3n10.dat
########## irp.py ##########

import localsolver
import sys
import math


def read_elem(filename):
    with open(filename) as f:
        return [str(elem) for elem in f.read().split()]


def main(instance_file, str_time_limit, sol_file):
    #
    # Reads instance data
    #
    (nb_customers, horizon_length, capacity, start_level_supplier, production_rate_supplier, holding_cost_supplier,
     start_level, max_level, demand_rate, holding_cost, distance_matrix, distance_supplier) = read_input_irp(
        instance_file)

    with localsolver.LocalSolver() as ls:
        #
        # Declares the optimization model
        #
        model = ls.model

        # Quantity of product delivered at each discrete time instant of the planning time horizon to each customer
        delivery = [[model.float(0, capacity) for i in range(nb_customers)] for t in range(horizon_length)]

        # Sequence of customers visited at each discrete time instant of the planning time horizon
        route = [model.list(nb_customers) for t in range(horizon_length)]

        # Customers receive products only if they are visited
        is_delivered = [[model.contains(route[t], i) for i in range(nb_customers)] for t in range(horizon_length)]

        # Create distance as an array to be able to access it with an "at" operator
        distance_array = model.array()
        for i in range(nb_customers):
            distance_array.add_operand(model.array(distance_matrix[i]))

        distance_supplier_array = model.array(distance_supplier)

        cost_route = [None for _ in range(horizon_length)]

        for t in range(horizon_length):
            sequence = route[t]
            c = model.count(sequence)

            # Distance traveled at instant t
            dist_selector = model.lambda_function(lambda i: model.at(distance_array, sequence[i - 1], sequence[i]))
            cost_route[t] = model.iif(c > 0, distance_supplier_array[sequence[0]] + model.sum(model.range(1, c),
                                                                                              dist_selector) +
                                      distance_supplier_array[sequence[c - 1]], 0)

        # Stockout constraints at the supplier
        inventory_supplier = [None for _ in range(horizon_length + 1)]
        inventory_supplier[0] = start_level_supplier
        for t in range(1, horizon_length + 1):
            inventory_supplier[t] = inventory_supplier[t - 1] - model.sum(
                delivery[t - 1][i] for i in range(nb_customers)) + production_rate_supplier
            if t != horizon_length:
                model.constraint(inventory_supplier[t] >= model.sum(delivery[t][i] for i in range(nb_customers)))

        # Stockout constraints at the customers
        inventory = [[None for _ in range(horizon_length + 1)] for _ in range(nb_customers)]
        for i in range(nb_customers):
            inventory[i][0] = start_level[i]
            for t in range(1, horizon_length + 1):
                inventory[i][t] = inventory[i][t - 1] + delivery[t - 1][i] - demand_rate[i]
                model.constraint(inventory[i][t] >= 0)

        for t in range(horizon_length):
            # Capacity constraints
            model.constraint(model.sum((delivery[t][i]) for i in range(nb_customers)) <= capacity)

            # Maximum level constraints
            for i in range(nb_customers):
                model.constraint(delivery[t][i] <= max_level[i] - inventory[i][t])
                model.constraint(delivery[t][i] <= max_level[i] * is_delivered[t][i])

        # Total inventory cost at the supplier
        total_cost_inventory_supplier = holding_cost_supplier * model.sum(
            inventory_supplier[t] for t in range(horizon_length + 1))

        # Total inventory cost at customers
        total_cost_inventory = model.sum(
            model.sum(holding_cost[i] * inventory[i][t] for t in range(horizon_length + 1)) for i in
            range(nb_customers))

        # Total transportation cost
        total_cost_route = model.sum(cost_route[t] for t in range(horizon_length))

        # Objective: minimize the sum of all costs
        objective = total_cost_inventory_supplier + total_cost_inventory + total_cost_route
        model.minimize(objective)

        model.close()

        #
        # Parameterizes the solver
        #
        ls.param.time_limit = int(str_time_limit)

        ls.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)
        #
        if len(sys.argv) >= 3:
            with open(sol_file, 'w') as f:
                f.write("%d\n" % (total_cost_route.value))
                for t in range(horizon_length):
                    for customer in route[t].value:
                        f.write("%d " % (customer + 1))
                    f.write("\n")


# The input files follow the "Archetti" format.
def read_input_irp(filename):
    file_it = iter(read_elem(filename))

    nb_customers = int(next(file_it)) - 1
    horizon_length = int(next(file_it))
    capacity = int(next(file_it))

    x_coord = [None] * nb_customers
    y_coord = [None] * nb_customers
    start_level = [None] * nb_customers
    max_level = [None] * nb_customers
    min_level = [None] * nb_customers
    demand_rate = [None] * nb_customers
    holding_cost = [None] * nb_customers

    next(file_it)
    x_coord_supplier = float(next(file_it))
    y_coord_supplier = float(next(file_it))
    start_level_supplier = int(next(file_it))
    production_rate_supplier = int(next(file_it))
    holding_cost_supplier = float(next(file_it))
    for i in range(nb_customers):
        next(file_it)
        x_coord[i] = float(next(file_it))
        y_coord[i] = float(next(file_it))
        start_level[i] = int(next(file_it))
        max_level[i] = int(next(file_it))
        min_level[i] = int(next(file_it))
        demand_rate[i] = int(next(file_it))
        holding_cost[i] = float(next(file_it))

    # Compute distance matrix
    distance_matrix = compute_distance_matrix(x_coord, y_coord)
    distance_supplier = compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord)

    return (
    nb_customers, horizon_length, capacity, start_level_supplier, production_rate_supplier, holding_cost_supplier,
    start_level, max_level, demand_rate, holding_cost, distance_matrix, distance_supplier)


# Computes the distance matrix
def compute_distance_matrix(x_coord, y_coord):
    nb_customers = len(x_coord)
    distance_matrix = [[None for i in range(nb_customers)] for j in range(nb_customers)]
    for i in range(nb_customers):
        distance_matrix[i][i] = 0
        for j in range(nb_customers):
            dist = compute_dist(x_coord[i], x_coord[j], y_coord[i], y_coord[j])
            distance_matrix[i][j] = dist
            distance_matrix[j][i] = dist
    return distance_matrix


# Computes the distances to the supplier
def compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord):
    nb_customers = len(x_coord)
    distance_supplier = [None] * nb_customers
    for i in range(nb_customers):
        dist = compute_dist(x_coord_supplier, x_coord[i], y_coord_supplier, y_coord[i])
        distance_supplier[i] = dist
    return distance_supplier


def compute_dist(xi, xj, yi, yj):
    return round(math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2)))


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python irp.py input_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    sol_file = sys.argv[2] if len(sys.argv) > 2 else None
    str_time_limit = sys.argv[3] if len(sys.argv) > 3 else "20"

    main(instance_file, str_time_limit, sol_file)
Compilation / Execution (Windows)
cl /EHsc irp.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver110.lib
irp instances\abs3n10.dat
Compilation / Execution (Linux)
g++ irp.cpp -I/opt/localsolver_11_0/include -llocalsolver110 -lpthread -o irp
./irp instances/abs3n10.dat
/********** irp.cpp **********/

#include <iostream>
#include <fstream>
#include <vector>
#include <cstring>
#include <cmath>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

class Irp {
public:
    // LocalSolver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Horizon length
    int horizonLength;

    // Capacity
    int capacity;

    // Start level at the supplier
    lsint startLevelSupplier;

    // Production rate of the supplier
    lsint productionRateSupplier;

    // Holding costs of the supplier
    lsdouble holdingCostSupplier;

    // Start level of the customers
    vector<lsint> startLevel;

    // Max level of the customers
    vector<lsint> maxLevel;

    // Demand rate of the customers
    vector<lsint> demandRate;

    // Holding costs of the customers
    vector<lsdouble> holdingCost;

    // Distance matrix
    vector<vector<lsint>> distanceMatrix;

    // Distance to depot
    vector<lsint> distanceSupplier;

    // Decision variables
    vector<vector<LSExpression>> delivery;

    // Decision variables
    vector<LSExpression> route;

    // Are the customers receiving products
    vector<vector<LSExpression>> isDelivered;

    // Total inventory cost at the supplier
    LSExpression totalCostInventorySupplier;

    // Total inventory cost at customers
    LSExpression totalCostInventory;

    // Total transportation cost
    LSExpression totalCostRoute;

    // Objective
    LSExpression objective;

    // Constructor
    Irp() {
    }

    // Reads instance data.
    void readInstance(const string& fileName) {
        readInputIrp(fileName);
    }

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

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

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

        // Create distance as an array to be able to access it with an "at" operator
        LSExpression distanceArray = model.array();
        for (int i = 0; i < nbCustomers; i++) {
            distanceArray.addOperand(model.array(distanceMatrix[i].begin(), distanceMatrix[i].end()));
        }
        LSExpression distanceSupplierArray = model.array(distanceSupplier.begin(), distanceSupplier.end());
        
        isDelivered.resize(horizonLength);
        vector<LSExpression> costRoute(horizonLength);

        for (int t = 0; t < horizonLength; t++) {
            LSExpression sequence = route[t];
            LSExpression c = model.count(sequence);
            isDelivered[t].resize(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.createLambdaFunction([&](LSExpression i) { return model.at(distanceArray, sequence[i - 1], sequence[i]); }); 
            costRoute[t] = model.iif(c > 0, distanceSupplierArray[sequence[0]] + model.sum(model.range(1, c), distSelector) + distanceSupplierArray[sequence[c - 1]], 0);
        }
        
        // Stockout constraints at the supplier
        vector<LSExpression> inventorySupplier(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].begin(), delivery[t-1].end()) + productionRateSupplier;
            if (t != horizonLength){
                model.constraint(inventorySupplier[t] >= model.sum(delivery[t].begin(), delivery[t].end()));
            }
        }

        // Stockout constraints at the customers
        vector<vector<LSExpression>> inventory(nbCustomers);
        for (int i = 0; i < nbCustomers; i++){
            inventory[i].resize(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].begin(), delivery[t].end()) <= 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.begin(), inventorySupplier.end());

        // Total inventory cost at customers
        vector<LSExpression> costInventoryCustomer(nbCustomers);
        for (int i = 0; i < nbCustomers; i++){
            costInventoryCustomer[i] = holdingCost[i] * model.sum(inventory[i].begin(), inventory[i].end());
        }
        totalCostInventory = model.sum(costInventoryCustomer.begin(), costInventoryCustomer.end());

        // Total transportation cost
        totalCostRoute = model.sum(costRoute.begin(), costRoute.end());
        
        // Objective: minimize the sum of all costs
        objective = 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(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        
        outfile << totalCostRoute.getValue() << endl;
        for (int t = 0; t < horizonLength; t++) {
            LSCollection routeCollection = route[t].getCollectionValue();
            for (lsint i = 0; i < routeCollection.count(); i++) {
                outfile << routeCollection[i] + 1 << " ";
            }
            outfile << endl;
        }        
    }

private:
    // The input files follow the "Archetti" format.
    void readInputIrp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }

        infile >> nbCustomers;
        nbCustomers -= 1;
        infile >> horizonLength;
        infile >> capacity;
        int idSupplier;
        double xCoordSupplier, yCoordSupplier;
        vector<int> id;
        vector<double> xCoord, yCoord;
        infile >> idSupplier;
        infile >> xCoordSupplier;
        infile >> yCoordSupplier;
        infile >> startLevelSupplier;
        infile >> productionRateSupplier;
        infile >> holdingCostSupplier;
        vector<int> minLevel;
        id.resize(nbCustomers);
        xCoord.resize(nbCustomers);
        yCoord.resize(nbCustomers);
        startLevel.resize(nbCustomers);
        maxLevel.resize(nbCustomers);
        minLevel.resize(nbCustomers);
        demandRate.resize(nbCustomers);
        holdingCost.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            infile >> id[i];
            infile >> xCoord[i];
            infile >> yCoord[i];
            infile >> startLevel[i];
            infile >> maxLevel[i];
            infile >> minLevel[i];
            infile >> demandRate[i];
            infile >> holdingCost[i];
        }
        
        printf("%lf", holdingCost[3]);

        // Compute distance matrix
        computeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);

        infile.close();
    }

    // Computes the distance matrix
    void computeDistanceMatrix(double xCoordSupplier, double yCoordSupplier, const vector<double>& xCoord, const vector<double>& yCoord) {
        distanceMatrix.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; i++) {
            distanceMatrix[i].resize(nbCustomers);
        }
        for (int i = 0; i < nbCustomers; i++) {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; j++) {
                lsint distance = computeDist(xCoord[i], xCoord[j], yCoord[i], yCoord[j]);
                distanceMatrix[i][j] = distance;
                distanceMatrix[j][i] = distance;
            }
        }
        
        distanceSupplier.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distanceSupplier[i] = computeDist(xCoordSupplier, xCoord[i], yCoordSupplier, yCoord[i]);
        }
    }

    lsint computeDist(double xi, double xj, double yi, double yj) {
        return round(sqrt(pow(xi - xj, 2) + pow(yi - yj, 2)));
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: irp inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }

    const char* instanceFile = argv[1];
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "20";

    try {
        Irp model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL) model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc Irp.cs /reference:localsolvernet.dll
Irp instances\abs3n10.dat
/********** 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))));
    }
}
Compilation / Execution (Windows)
javac Irp.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Irp instances/abs3n10.dat
Compilation / Execution (Linux)
javac Irp.java -cp /opt/localsolver_11_0/bin/localsolver.jar
java -cp /opt/localsolver_11_0/bin/localsolver.jar:. Irp instances/abs3n10.dat
/********** Irp.java **********/

import java.util.*;
import java.io.*;
import localsolver.*;

public class Irp {
    // Solver
    private final LocalSolver localsolver;

    // Number of customers
    private int nbCustomers;

    // Horizon length
    private int horizonLength;

    // Capacity
    private int capacity;

    // Start level at the supplier
    private int startLevelSupplier;

    // Production rate of the supplier
    private int productionRateSupplier;

    // Holding costs of the supplier
    private double holdingCostSupplier;

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

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

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

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

    // Distance matrix
    private long[][] distanceMatrix;
    
    // Distances between customers and supplier
    private 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
    private LSExpression[] costRoute;

    // Inventory at the supplier
    private LSExpression[] inventorySupplier;

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

    // Total inventory cost at the supplier
    private LSExpression totalCostInventorySupplier;

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

    // Total inventory cost at customers
    private LSExpression totalCostInventory;

    // Total transportation cost
    private LSExpression totalCostRoute;

    // Objective
    private LSExpression objective;

    private Irp(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    private void solve(int limit) {
        // Declares the optimization model.
        LSModel model = localsolver.getModel();
        
        delivery = new LSExpression[horizonLength][nbCustomers];
        route =  new LSExpression[horizonLength];
        isDelivered = new LSExpression[horizonLength][nbCustomers];
        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++){
            for (int i = 0; i < nbCustomers; i++){
                delivery[t][i] = model.floatVar(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.listVar(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);
            
            // 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 -> model.at(
                        distanceArray,
                        model.at(sequence, model.sub(i, 1)),
                        model.at(sequence, i)));
            costRoute[t] = model.iif(model.gt(c, 0), model.sum(model.sum(
                                model.at(distanceSupplierArray, model.at(sequence, 0)), model.sum(model.range(1, c), distSelector)), 
                                model.at(distanceSupplierArray, model.at(sequence, model.sub(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] = model.sum(model.sub(inventorySupplier[t-1], model.sum(delivery[t-1])), productionRateSupplier);
            if (t != horizonLength){
                model.constraint(model.geq(inventorySupplier[t], model.sum(delivery[t])));
            }
        }

        // Stockout constraints at the customers
        inventory = new LSExpression[nbCustomers][horizonLength+1];
        for (int i = 0; i < nbCustomers; i++){
            inventory[i][0] = model.createConstant(startLevel[i]);
            for (int t = 1; t < horizonLength + 1; t++){
                inventory[i][t] = model.sub(model.sum(inventory[i][t-1], delivery[t-1][i]), demandRate[i]);
                model.constraint(model.geq(inventory[i][t], 0));
            }
        }

        for (int t = 0; t < horizonLength; t++){
            // Capacity constraints
            model.constraint(model.leq(model.sum(delivery[t]), capacity));

            // Maximum level constraints
            for (int i = 0; i < nbCustomers; i++){
                model.constraint(model.leq(delivery[t][i], model.sub(maxLevel[i], inventory[i][t])));
                model.constraint(model.leq(delivery[t][i], model.prod(maxLevel[i], isDelivered[t][i])));
            }
        }

        // Total inventory cost at the supplier
        totalCostInventorySupplier = model.prod(holdingCostSupplier, model.sum(inventorySupplier));

        // Total inventory cost at customers
        costInventoryCustomer = new LSExpression[nbCustomers];
        for (int i = 0; i < nbCustomers; i++){
            costInventoryCustomer[i] = model.prod(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)
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(totalCostRoute.getValue());
            for (int t = 0; t < horizonLength; t++) {
                LSCollection routeCollection = route[t].getCollectionValue();
                for (int i = 0; i < routeCollection.count(); i++) {
                    output.print((routeCollection.get(i) + 1) + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "Archetti" format.
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            String[] splitted;
            splitted = input.nextLine().split("\\s+");
            nbCustomers = Integer.parseInt(splitted[0]) - 1;
            horizonLength = Integer.parseInt(splitted[1]);
            capacity = Integer.parseInt(splitted[2]);
            splitted = input.nextLine().split("\\s+");
            double xCoordSupplier = Double.parseDouble(splitted[2]);
            double yCoordSupplier = Double.parseDouble(splitted[3]);
            startLevelSupplier = Integer.parseInt(splitted[4]);
            productionRateSupplier = Integer.parseInt(splitted[5]);
            holdingCostSupplier = Double.parseDouble(splitted[6]);
            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.nextLine().split("\\s+");
                xCoord[i] = Double.parseDouble(splitted[2]);
                yCoord[i] = Double.parseDouble(splitted[3]);
                startLevel[i] = Integer.parseInt(splitted[4]);
                maxLevel[i] = Integer.parseInt(splitted[5]);
                demandRate[i] = Integer.parseInt(splitted[7]);
                holdingCost[i] = Double.parseDouble(splitted[8]);
            }
            // 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][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 Math.round(Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2)));
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java Irp inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        try (LocalSolver localsolver = new LocalSolver()) {
            String instanceFile = args[0];
            String outputFile = args.length > 1 ? args[1] : null;
            String strTimeLimit = args.length > 2 ? args[2] : "20";

            Irp model = new Irp(localsolver);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch(Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}