# 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¶

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.

## 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;

// 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] +
sum(1..c-1, i => distance[sequence[i-1]][sequence[i]]) +
distanceSupplier[sequence[c-1]] : 0;
}

// Stockout constraints at the supplier
inventorySupplier <- 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] <- 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();
}
}

// General data
nbCustomers = firstLine.toInt()-1;
horizonLength = firstLine.toInt();
capacity = firstLine.toInt();

// Supplier data
xCoordSupplier = secondLine.toDouble();
yCoordSupplier = secondLine.toDouble();
startLevelSupplier = secondLine.toInt();
productionRate_supplier = secondLine.toInt();
holdingCostSupplier = secondLine.toDouble();

// Customers data
for[i in 0..nbCustomers-1] {
xCoord[i] = elements.toDouble();
yCoord[i] = elements.toDouble();
startLevel[i] = elements.toInt();
maxLevel[i] = elements.toInt();
minLevel[i] = elements.toInt();
demandRate[i] = elements.toInt();
holdingCost[i] = elements.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

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

def main(instance_file, str_time_limit, sol_file):
#
#
(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_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] + 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 = 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] = 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.

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
sol_file = sys.argv if len(sys.argv) > 2 else None
str_time_limit = sys.argv 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() {
}

}

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++) {
}
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] + model.sum(model.range(1, c), distSelector) + distanceSupplierArray[sequence[c - 1]], 0);
}

// Stockout constraints at the supplier
vector<LSExpression> inventorySupplier(horizonLength+1);
inventorySupplier = 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] = 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.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.
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);

// 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;
const char* solFile = argc > 2 ? argv : NULL;
const char* strTimeLimit = argc > 3 ? argv : "20";

try {
Irp model;
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] + model.Sum(model.Range(1, c), distSelector) + distanceSupplierArray[sequence[c - 1]], 0);
}

// Stockout constraints at the supplier
inventorySupplier =  new LSExpression[horizonLength+1];
inventorySupplier = 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] = 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;
string outputFile = args.Length > 1 ? args : null;
string strTimeLimit = args.Length > 2 ? args : "20";

using (Irp model = new Irp())
{
model.Solve(int.Parse(strTimeLimit));
if (outputFile != null)
model.WriteSolution(outputFile);
}
}

// The input files follow the "Archetti" format.
{
{
string[] splitted;
splitted = input.ReadLine().Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
nbCustomers = int.Parse(splitted) - 1;
horizonLength = int.Parse(splitted);
capacity = int.Parse(splitted);
splitted = input.ReadLine().Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
Console.WriteLine(1);
Console.WriteLine(splitted);
double xCoordSupplier = double.Parse(splitted, CultureInfo.InvariantCulture);
Console.WriteLine(2);
double yCoordSupplier = double.Parse(splitted, CultureInfo.InvariantCulture);
Console.WriteLine(3);
startLevelSupplier = int.Parse(splitted);
Console.WriteLine(4);
productionRateSupplier = int.Parse(splitted);
holdingCostSupplier = double.Parse(splitted, 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, CultureInfo.InvariantCulture);
yCoord[i] = double.Parse(splitted, CultureInfo.InvariantCulture);
startLevel[i] = int.Parse(splitted);
maxLevel[i] = int.Parse(splitted);
demandRate[i] = int.Parse(splitted);
holdingCost[i] = double.Parse(splitted, 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 = 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] = 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) - 1;
horizonLength = Integer.parseInt(splitted);
capacity = Integer.parseInt(splitted);
splitted = input.nextLine().split("\\s+");
double xCoordSupplier = Double.parseDouble(splitted);
double yCoordSupplier = Double.parseDouble(splitted);
startLevelSupplier = Integer.parseInt(splitted);
productionRateSupplier = Integer.parseInt(splitted);
holdingCostSupplier = Double.parseDouble(splitted);
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);
yCoord[i] = Double.parseDouble(splitted);
startLevel[i] = Integer.parseInt(splitted);
maxLevel[i] = Integer.parseInt(splitted);
demandRate[i] = Integer.parseInt(splitted);
holdingCost[i] = Double.parseDouble(splitted);
}
// 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;
String outputFile = args.length > 1 ? args : null;
String strTimeLimit = args.length > 2 ? args : "20";

Irp model = new Irp(localsolver);