# Split Delivery Vehicle Routing (SDVRP)¶

## Principles learned¶

• Add multiple list decision variables

• Combine list and float decision variables

• Define a sequence of expressions

## Problem¶

In the split delivery vehicle routing problem (SDVRP), a fleet of delivery vehicles with uniform capacity must service customers with known demand for a single commodity. The vehicles start and end their routes at a common depot. Each customer can be served by multiple vehicles. A sequence of customers is assigned to each truck of the fleet. The objective is to minimize the total distance traveled by the trucks while all customers are served at least their demand and the total demand served by each truck does not exceed its own capacity.

## Data¶

The instances provided were introduced by Belenguer et at. instances and use the format of the DIMACS challenge.

The format of the data files is as follows:

• The first line contains the number of customers and the capacity Q of each truck

• The second line contains the demand of the customers

• The rest of the lines represents the coordinates of the warehouse followed by the coordinates of the customers

## Program¶

This LocalSolver model is a relaxation of the CVRP model. We refer the reader to this model for the routing aspect of the problem. In the split delivery routing problem, it’s no longer required to deliver customer demand by a single truck, meaning that the demand can be splitted and served by different trucks. Dror and Trudeau (1989) . showed that important cost savings can be obtained using split deliveries compared to customer single delivery.

While in the CVRP model the “partition” operator operator is used to ensure that each customer is visited exactly once, here we use the “cover” operator to ensure that all customers must be visited at least once, without any constraint on the number of trucks that visits each customer.

The quantity delievered by the truck `k` to the customer `i` if visited is `quantity[k][i]`. The total quantity delivered by each truck is computed thanks to the `sum` operator over all customers visited by a truck. This quantity is constrained to be lower than the truck capacity.

The quantity delivered to a customer is computed by summing the quantity delivered by all the trucks visiting this customer. This quantity is constrained to satisfy the customer demand.

Execution:
localsolver sdvrp.lsp inFileName=instances/S51D4.sd [lsTimeLimit=] [solFileName=]
```use io;

function input() {
usage = "Usage: localsolver sdvrp.lsp " +
"inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]";

if (inFileName == nil) throw usage;

// Compute distance matrix
computeDistanceMatrix();
}

/* Declare the optimization model */
function model () {
// Sequence of customers visited by each truck

// All customers must be visited at least by one truck

// Quantity carried by each truck for each customer
quantity[k in 0...nbTrucks][i in 0...nbCustomers] <- float(0, demands[i]);

for [k in 0...nbTrucks] {
local c <- count(sequence);

// The quantity needed in each route must not exceed the truck capacity
routeQuantity[k] <- sum(sequence, j => quantity[k][j]);
constraint routeQuantity[k] <= truckCapacity;

// A truck is used if it visits at least one customer
trucksUsed[k] <- c > 0;

// Distance traveled by truck k
routeDistances[k] <- sum(1...c,
i => distanceMatrix[sequence[i-1]][sequence[i]]) + (trucksUsed[k] ?
(distanceDepot[sequence] + distanceDepot[sequence[c-1]]) :
0);
}

for [i in 0...nbCustomers] {
// Each customer must receive at least its demand
quantityServed[i] <- sum[k in 0...nbTrucks](quantity[k][i]
constraint quantityServed[i] >= demands[i];
}

totalDistance <- sum[k in 0...nbTrucks](routeDistances[k]);

// Objective: minimize the distance traveled
minimize totalDistance;
}

/* Write the solution in a file with the following format:
* - Each line k contains the customers visited by the truck k */
function output () {
if (solFileName == nil) return;
local outfile = io.openWrite(solFileName);

outfile.println(totalDistance.value);
for [customer in 0...nbCustomers] {
outfile.print(customer + 1, ": ");
for [k in 0...nbTrucks] {
if (trucksUsed[k].value != 1) continue;
outfile.print(k + 1, " ");
}
outfile.println();
}
}

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

nbTrucks = nbCustomers;

// Extracting the demand of each customer

// Extracting the coordinates of the depot and the customers
for [i in 0...nbCustomers + 1] {
}
}

function computeDistanceMatrix() {
// Computing the distance between customers
for [i in 0...nbCustomers] {
distanceMatrix[i][i] = 0;
for [j in i+1...nbCustomers] {
// +1 because computeDist expected original data indices,
// with customers in 1...nbCustomers+1
local localDistance = computeDist(i+1, j+1);
distanceMatrix[j][i] = localDistance;
distanceMatrix[i][j] = localDistance;
}
}
// Computing the distance between depot and customers
for [i in 0...nbCustomers] {
local localDistance = computeDist(0, i+1);
distanceDepot[i] = localDistance;
}
}

function computeDist(i, j) {
local exactDist = sqrt(pow((nodesX[i] - nodesX[j]), 2) + pow((nodesY[i] - nodesY[j]), 2));
return floor(exactDist + 0.5);
}
```
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python sdvrp.py instances\S51D4.sd
Execution (Linux)
export PYTHONPATH=/opt/localsolver_12_0/bin/python
python sdvrp.py instances/S51D4.sd
```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, truck_capacity, dist_matrix_data, dist_depot_data, demands_data = \
nb_trucks = nb_customers

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

# Sequence of customers visited by each truck
customers_sequences = [model.list(nb_customers) for _ in range(nb_trucks)]

# Quantity carried by each truck for each customer
quantity = [None] * nb_trucks
for k in range(nb_trucks):
quantity[k] = [model.float(0, demands_data[i]) for i in range(nb_customers)]

# All customers must be visited at least by one truck
model.constraint(model.cover(customers_sequences))

# Create LocalSolver arrays to be able to access them with "at" operators
dist_matrix = model.array(dist_matrix_data)
dist_depot = model.array(dist_depot_data)

route_distances = [None] * nb_trucks
trucks_used = [None] * nb_trucks

for k in range(nb_trucks):
sequence = customers_sequences[k]
c = model.count(sequence)

# A truck is used if it visits at least one customer
trucks_used[k] = model.count(sequence) > 0

# The quantity carried in each route must not exceed the truck capacity
quantity_array = model.array(quantity[k])
quantity_lambda = model.lambda_function(lambda j: quantity_array[j])
route_quantity = model.sum(sequence, quantity_lambda)
model.constraint(route_quantity <= truck_capacity)

# Distance traveled by each truck
dist_lambda = model.lambda_function(
lambda i: model.at(dist_matrix, sequence[i - 1], sequence[i]))
route_distances[k] = model.sum(model.range(1, c), dist_lambda) \
+ model.iif(
trucks_used[k],
dist_depot[sequence] + dist_depot[sequence[c - 1]],
0)

for i in range(nb_customers):
# Each customer must receive at least its demand
quantity_served = model.sum(
quantity[k][i] * model.contains(customers_sequences[k], i)
for k in range(nb_trucks))
model.constraint(quantity_served >= demands_data[i])

total_distance = model.sum(route_distances)

# Objective: minimize the distance traveled
model.minimize(total_distance)

model.close()

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

ls.solve()

#
# Write the solution in a file with the following format:
# each line k contain the customers visited by the truck k
#
if len(sys.argv) >= 3:
with open(sol_file, 'w') as f:
f.write("%d \n" % total_distance.value)
for k in range(nb_trucks):
if trucks_used[k].value != 1:
continue
# Values in sequence are in 0...nbCustomers. +1 is to put it back in 1...nbCustomers+1
for customer in customers_sequences[k].value:
f.write("%d " % (customer + 1))
f.write("\n")

nb_customers = int(next(file_it))
capacity = int(next(file_it))

demands = [None] * nb_customers
for i in range(nb_customers):
demands[i] = int(next(file_it))

# Extracting the coordinates of the depot and the customers
customers_x = [None] * nb_customers
customers_y = [None] * nb_customers
depot_x = float(next(file_it))
depot_y = float(next(file_it))
for i in range(nb_customers):
customers_x[i] = float(next(file_it))
customers_y[i] = float(next(file_it))

distance_matrix = compute_distance_matrix(customers_x, customers_y)
distance_depots = compute_distance_depots(depot_x, depot_y, customers_x, customers_y)
return nb_customers, capacity, distance_matrix, distance_depots, demands

# Compute the distance between two customers
def compute_distance_matrix(customers_x, customers_y):
nb_customers = len(customers_x)
distance_matrix = [[None for _ in range(nb_customers)] for _ in range(nb_customers)]
for i in range(nb_customers):
distance_matrix[i][i] = 0
for j in range(nb_customers):
dist = compute_dist(customers_x[i], customers_x[j], customers_y[i], customers_y[j])
distance_matrix[i][j] = dist
distance_matrix[j][i] = dist
return distance_matrix

# Compute the distances to depot
def compute_distance_depots(depot_x, depot_y, customers_x, customers_y):
nb_customers = len(customers_x)
distance_depots = [None] * nb_customers
for i in range(nb_customers):
dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
distance_depots[i] = dist
return distance_depots

def compute_dist(xi, xj, yi, yj):
exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
return int(math.floor(exact_dist + 0.5))

if __name__ == '__main__':
if len(sys.argv) < 2:
print("Usage: python sdvrp.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 sdvrp.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver120.lib
sdvrp instances\S51D4.sd
Compilation / Execution (Linux)
g++ sdvrp.cpp -I/opt/localsolver_12_0/include -llocalsolver120 -lpthread -o sdvrp
./sdvrp instances/S51D4.sd
```#include "localsolver.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace localsolver;
using namespace std;

class Sdvrp {
public:
// LocalSolver
LocalSolver localsolver;

// Number of customers
int nbCustomers;

// Capacity of the trucks
int truckCapacity;

// Demand of each customer
vector<int> demandsData;

// Distance matrix
vector<vector<int>> distMatrixData;

// Distance between customers and depot
vector<int> distDepotData;

// Number of trucks
int nbTrucks;

// Decision variables

// Are the trucks actually used
vector<LSExpression> trucksUsed;

// Distance traveled by all the trucks
LSExpression totalDistance;

// Quantity carried by each truck
vector<vector<LSExpression>> quantity;

nbTrucks = nbCustomers;
}

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

// Sequence of customers visited by each truck
for (int k = 0; k < nbTrucks; ++k) {
}

// Quantity carried by each truck for each customer
quantity.resize(nbTrucks);
for (int k = 0; k < nbTrucks; ++k) {
quantity[k].resize(nbTrucks);
for (int i = 0; i < nbCustomers; ++i) {
quantity[k][i] = model.floatVar(0, demandsData[i]);
}
}

// All customers must be visited by exactly one truck

// Create LocalSolver arrays to be able to access them with "at" operators
LSExpression distMatrix = model.array();
for (int n = 0; n < nbCustomers; ++n) {
}
LSExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());

trucksUsed.resize(nbTrucks);
vector<LSExpression> distRoutes(nbTrucks);

for (int k = 0; k < nbTrucks; ++k) {
LSExpression c = model.count(sequence);

// A truck is used if it visits at least one customer
trucksUsed[k] = c > 0;

// The quantity carried in each route must not exceed the truck capacity
LSExpression quantityArray = model.array(quantity[k].begin(), quantity[k].end());
LSExpression quantityLambda = model.createLambdaFunction(
[&](LSExpression j) { return quantityArray[j]; });
LSExpression routeQuantity = model.sum(sequence, quantityLambda);
model.constraint(routeQuantity <= truckCapacity);

// Distance traveled by truck k
LSExpression distLambda = model.createLambdaFunction([&](LSExpression i) {
return model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i));
});
distRoutes[k] = model.sum(model.sum(model.range(1, c), distLambda),
model.iif(trucksUsed[k],
model.sum(model.at(distDepot, model.at(sequence, 0)),
model.at(distDepot, model.at(sequence, c - 1))),
0));
}

for (int i = 0; i < nbCustomers; ++i) {
// Each customer must receive at least its demand
LSExpression quantityServed = model.sum();
for (int k = 0; k < nbTrucks; ++k)
model.constraint(quantityServed >= demandsData[i]);
}

totalDistance = model.sum(distRoutes.begin(), distRoutes.end());

// Objective: minimize the distance traveled
model.minimize(totalDistance);

model.close();

// Parameterize the solver
localsolver.getParam().setTimeLimit(limit);

localsolver.solve();
}

/* Write the solution in a file with the following format:
*  - Each line k contains the customers visited by the truck k */
void writeSolution(const string& fileName) {
ofstream outfile;
outfile.open(fileName.c_str());

outfile << totalDistance.getValue() << endl;
for (int k = 0; k < nbTrucks; ++k) {
if (trucksUsed[k].getValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +1 is to put it back in 1...nbCustomers+1
for (int i = 0; i < customersCollection.count(); ++i) {
outfile << customersCollection[i] + 1 << " ";
}
outfile << endl;
}
}

private:
ifstream infile(fileName.c_str());
if (!infile.is_open()) {
throw std::runtime_error("File cannot be opened.");
}

string str;
char* line;
istringstream iss;
getline(infile, str);
line = strdup(str.c_str());
iss.str(line);
iss >> nbCustomers;
iss >> truckCapacity;
nbTrucks = nbCustomers;
getline(infile, str);
line = strdup(str.c_str());
istringstream stream(line);
demandsData.resize(nbCustomers);
int demand;
for (int n = 0; n < nbCustomers; ++n) {
stream >> demand;
demandsData[n] = demand;
}
vector<int> customersX(nbCustomers);
vector<int> customersY(nbCustomers);
int depotX, depotY;
getline(infile, str);
line = strdup(str.c_str());
istringstream iss_1(line);
iss_1 >> depotX;
iss_1 >> depotY;
for (int i = 0; i < nbCustomers; ++i) {
getline(infile, str);
line = strdup(str.c_str());
istringstream iss(line);
iss >> customersX[i];
iss >> customersY[i];
}

computeDistanceMatrix(depotX, depotY, customersX, customersY);

infile.close();
}

// Compute the distance matrix
void computeDistanceMatrix(int depotX, int depotY, const vector<int>& customersX, const vector<int>& customersY) {
distMatrixData.resize(nbCustomers);
for (int i = 0; i < nbCustomers; ++i) {
distMatrixData[i].resize(nbCustomers);
}
for (int i = 0; i < nbCustomers; ++i) {
distMatrixData[i][i] = 0;
for (int j = i + 1; j < nbCustomers; ++j) {
int distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
distMatrixData[i][j] = distance;
distMatrixData[j][i] = distance;
}
}

distDepotData.resize(nbCustomers);
for (int i = 0; i < nbCustomers; ++i) {
distDepotData[i] = computeDist(depotX, customersX[i], depotY, customersY[i]);
}
}

static int computeDist(int xi, int xj, int yi, int yj) {
double exactDist = sqrt(pow((double)xi - xj, 2) + pow((double)yi - yj, 2));
return floor(exactDist + 0.5);
}
};

int main(int argc, char** argv) {
if (argc < 2) {
cerr << "Usage: Sdvrp 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 {
Sdvrp 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 Sdvrp.cs /reference:localsolvernet.dll
Sdvrp instances\S51D4.sd
```/********** Sdvrp.cs **********/
using System;
using System.IO;
using localsolver;

public class Sdvrp : IDisposable
{
// LocalSolver
LocalSolver localsolver;

// Number of customers
int nbCustomers;

// Capacity of the trucks
int truckCapacity;

// Demand of each customer
long[] demandsData;

// Distance matrix
long[][] distMatrixData;

// Distances between customers and depot
long[] distDepotData;

// Number of trucks
int nbTrucks;

// Decision variables

// Are the trucks actually used
LSExpression[] trucksUsed;

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

// Distance traveled by all the trucks
LSExpression totalDistance;

// Quantity carried by each truck
LSExpression[][] quantity;

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

{
}

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

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

trucksUsed = new LSExpression[nbTrucks];
distRoutes = new LSExpression[nbTrucks];
quantity = new LSExpression[nbTrucks][];

// Sequence of customers visited by each truck
for (int k = 0; k < nbTrucks; ++k)

// Quantity carried by each truck for each customer
for (int k = 0; k < nbTrucks; ++k)
{
quantity[k] = new LSExpression[nbCustomers];
for (int i = 0; i < nbCustomers; ++i)
quantity[k][i] = model.Float(0, demandsData[i]);
}

// Every customer must be visited by at least one truck

// Create LocalSolver arrays to be able to access them with "at" operators
LSExpression distMatrix = model.Array(distMatrixData);
LSExpression distDepot = model.Array(distDepotData);

for (int k = 0; k < nbTrucks; ++k)
{
LSExpression c = model.Count(sequence);

// A truck is used if it visits at least one customer
trucksUsed[k] = c > 0;

// The quantity carried in each route must not exceed the truck capacity
LSExpression quantityArray = model.Array(quantity[k]);
LSExpression quantityLambda = model.LambdaFunction(j => quantityArray[j]);
LSExpression routeQuantity = model.Sum(sequence, quantityLambda);
model.Constraint(routeQuantity <= truckCapacity);

// Distance traveled by truck k
LSExpression distLambda = model.LambdaFunction(
i => distMatrix[sequence[i - 1], sequence[i]]
);
distRoutes[k] =
model.Sum(model.Range(1, c), distLambda)
+ model.If(trucksUsed[k], distDepot[sequence] + distDepot[sequence[c - 1]], 0);
}

for (int i = 0; i < nbCustomers; ++i)
{
// Each customer must receive at least its demand
LSExpression quantityServed = model.Sum();
for (int k = 0; k < nbTrucks; ++k)
);
model.Constraint(quantityServed >= demandsData[i]);
}

totalDistance = model.Sum(distRoutes);

// Objective: minimize the distance traveled
model.Minimize(totalDistance);

model.Close();

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

localsolver.Solve();
}

/* Write the solution in a file with the following format:
*  - Each line k contains the customers visited by the truck k */
void WriteSolution(string fileName)
{
using (StreamWriter output = new StreamWriter(fileName))
{
output.WriteLine(totalDistance.GetValue());
for (int k = 0; k < nbTrucks; ++k)
{
if (trucksUsed[k].GetValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +1 is to put it back in 1...nbCustomers+1
for (int i = 0; i < customersCollection.Count(); ++i)
output.Write((customersCollection[i] + 1) + " ");
output.WriteLine();
}
}
}

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

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

{
if (line == null)
return new string;
return line.Split(new[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
}

{
{
string[] splitted;
splitted = SplitInput(input);

nbCustomers = int.Parse(splitted);
nbTrucks = nbCustomers;
truckCapacity = int.Parse(splitted);

demandsData = new long[nbCustomers];

splitted = SplitInput(input);
for (int n = 0; n < nbCustomers; ++n)
demandsData[n] = int.Parse(splitted[n]);

splitted = SplitInput(input);
int depotX = int.Parse(splitted);
int depotY = int.Parse(splitted);

int[] customersX = new int[nbCustomers];
int[] customersY = new int[nbCustomers];
for (int n = 0; n < nbCustomers; ++n)
{
splitted = SplitInput(input);
customersX[n] = int.Parse(splitted);
customersY[n] = int.Parse(splitted);
}

ComputeDistanceMatrix(depotX, depotY, customersX, customersY);
}
}

// Compute the distance matrix
private void ComputeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY)
{
distMatrixData = new long[nbCustomers][];
for (int i = 0; i < nbCustomers; ++i)
distMatrixData[i] = new long[nbCustomers];

for (int i = 0; i < nbCustomers; ++i)
{
distMatrixData[i][i] = 0;
for (int j = i + 1; j < nbCustomers; ++j)
{
long dist = ComputeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
distMatrixData[i][j] = dist;
distMatrixData[j][i] = dist;
}
}
distDepotData = new long[nbCustomers];
for (int i = 0; i < nbCustomers; ++i)
distDepotData[i] = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
}

private static long ComputeDist(int xi, int xj, int yi, int yj)
{
double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
return (long)Math.Floor(exactDist + 0.5);
}
}
```
Compilation / Execution (Windows)
javac Sdvrp.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Sdvrp instances\S51D4.sd
Compilation / Execution (Linux)
javac Sdvrp.java -cp /opt/localsolver_12_0/bin/localsolver.jar
java -cp /opt/localsolver_12_0/bin/localsolver.jar:. Sdvrp instances/S51D4.sd
```import java.util.*;
import java.io.*;
import localsolver.*;

public class Sdvrp {
// LocalSolver
private final LocalSolver localsolver;

// Number of customers
int nbCustomers;

// Capacity of the trucks
private int truckCapacity;

// Demand of each customer
private long[] demandsData;

// Distance matrix
private long[][] distMatrixData;

// Distances between customers and depot
private long[] distDepotData;

// Number of trucks
private int nbTrucks;

// Decision variables

// Are the trucks actually used
private LSExpression[] trucksUsed;

// Quantity carried by each truck
private LSExpression[][] quantity;

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

// Distance traveled by all the trucks
private LSExpression totalDistance;

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

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

trucksUsed = new LSExpression[nbTrucks];
distRoutes = new LSExpression[nbTrucks];
quantity = new LSExpression[nbTrucks][nbCustomers];

// Sequence of customers visited by each truck
for (int k = 0; k < nbTrucks; ++k)

// Quantity carried by each truck for each customer
for (int k = 0; k < nbTrucks; ++k)
for (int i = 0; i < nbCustomers; ++i)
quantity[k][i] = model.floatVar(0, demandsData[i]);

// Every customer must be visited by at least one truck

// Create LocalSolver arrays to be able to access them with "at" operators
LSExpression distMatrix = model.array(distMatrixData);
LSExpression distDepot = model.array(distDepotData);

for (int k = 0; k < nbTrucks; ++k) {
LSExpression c = model.count(sequence);

// A truck is used if it visits at least one customer
trucksUsed[k] = model.gt(c, 0);

// The quantity carried in each route must not exceed the truck capacity
LSExpression quantityArray = model.array(quantity[k]);
LSExpression quantityLambda = model.lambdaFunction(j -> model.at(quantityArray, j));
LSExpression routeQuantity = model.sum(sequence, quantityLambda);
model.constraint(model.leq(routeQuantity, truckCapacity));

// Distance traveled by truck k
LSExpression distLambda = model
.lambdaFunction(i -> model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)));
distRoutes[k] = model.sum(model.sum(model.range(1, c), distLambda),
model.iif(trucksUsed[k], model.sum(model.at(distDepot, model.at(sequence, 0)),
model.at(distDepot, model.at(sequence, model.sub(c, 1)))), 0));
}

for (int i = 0; i < nbCustomers; ++i) {
// Each customer must receive at least its demand
LSExpression quantityServed = model.sum();
for (int k = 0; k < nbTrucks; ++k)
model.constraint(model.geq(quantityServed, demandsData[i]));
}

totalDistance = model.sum(distRoutes);

// Objective: minimize the distance traveled
model.minimize(totalDistance);

model.close();

// Parameterize the solver
localsolver.getParam().setTimeLimit(limit);

localsolver.solve();
}

private void readInstance(String fileName) throws IOException {
try (Scanner input = new Scanner(new File(fileName))) {
nbCustomers = input.nextInt();
truckCapacity = input.nextInt();
nbTrucks = nbCustomers;
input.nextLine();
demandsData = new long[nbCustomers];
for (int n = 1; n <= nbCustomers; ++n) {
demandsData[n - 1] = input.nextInt();
}

input.nextLine();

int[] customersX = new int[nbCustomers];
int[] customersY = new int[nbCustomers];
int depotX = input.nextInt();
int depotY = input.nextInt();
input.nextLine();
for (int n = 1; n <= nbCustomers; ++n) {
customersX[n - 1] = input.nextInt();
customersY[n - 1] = input.nextInt();
input.nextLine();
}

computeDistanceMatrix(depotX, depotY, customersX, customersY);
}
}

// Compute the distance matrix
private void computeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY) {
distMatrixData = new long[nbCustomers][nbCustomers];
for (int i = 0; i < nbCustomers; ++i) {
distMatrixData[i][i] = 0;
for (int j = i + 1; j < nbCustomers; ++j) {
long dist = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
distMatrixData[i][j] = dist;
distMatrixData[j][i] = dist;
}
}

distDepotData = new long[nbCustomers];
for (int i = 0; i < nbCustomers; ++i) {
distDepotData[i] = computeDist(depotX, customersX[i], depotY, customersY[i]);
}
}

private static long computeDist(int xi, int xj, int yi, int yj) {
double exactDist = Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
return (long) Math.floor(exactDist + 0.5);
}

public static void main(String[] args) {
if (args.length < 1) {
System.err.println("Usage: java Sdvrp 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";

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

/* Write the solution in a file with the following format:
* - Each line k contains the customers visited by the truck k */
private void writeSolution(String fileName) throws IOException {
try (PrintWriter output = new PrintWriter(fileName)) {
output.println(totalDistance.getValue());
for (int k = 0; k < nbTrucks; ++k) {
if (trucksUsed[k].getValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +1 is to put it back in 1...nbCustomers+1