# Split Delivery Vehicle Routing (SDVRP)¶

## Principles learned¶

• Add multiple list decision variables

• Add a cover constraint

• 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=]
```/********** sdvrp.lsp **********/
use io;

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

if (inFileName == nil) throw usage;

// Compute distance matrix
computeDistanceMatrix();
}

/* Declares the optimization model */
function model () {
// Sequence of customers visited by each truck
customersSequences[k in 0..nbTrucks-1] <- list(nbCustomers);

// All customers must be visited at least by one truck
constraint cover[k in 0..nbTrucks-1](customersSequences[k]);

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

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

// The quantity needed in each route must not exceed the truck capacity
routeQuantity[k] <- sum[i in 0..nbCustomers-1](quantity[k][i] * contains(sequence, i));
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-1,
i => distanceMatrix[sequence[i-1]][sequence[i]]) + (trucksUsed[k] ?
(distanceWarehouse[sequence] + distanceWarehouse[sequence[c-1]]) :
0);
}

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

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

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

function output () {
if (solFileName == nil) return;
local outfile = io.openWrite(solFileName);

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

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

local inFile = io.openRead(inFileName);
nbTrucks = nbCustomers;

// Extracting the demand of each customer
demands[i in 0..nbCustomers-1] = inFile.readInt();

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

/* Computes the distance between customers */
function computeDistanceMatrix() {
// Computing the distance between customers
for[i in 0..nbCustomers-1] {
distanceMatrix[i][i] = 0;
for[j in i+1..nbCustomers-1] {
// +1 because computeDist expected original data indices,
// with customers in 1..nbCustomers
local localDistance = computeDist(i+1, j+1);
distanceMatrix[j][i] = localDistance;
distanceMatrix[i][j] = localDistance;
}
}
// Computing the distance between warehouse and customers
for[i in 0..nbCustomers-1] {
local localDistance = computeDist(0, i+1);
distanceWarehouse[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_11_0/bin/python
python sdvrp.py instances/S51D4.sd
```########## sdvrp.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):
#
# Reads instance data
#
(nb_customers, truck_capacity, distance_matrix, distance_warehouses, demands) = read_input_sdvrp(instance_file)
nb_trucks = nb_customers

with localsolver.LocalSolver() as ls:
#
# Declares 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[i]) for i in range(nb_customers)]

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

# Creates distance as an array to be able to access it with an "at" operator
distance_array = model.array(distance_matrix)
distance_warehouse_array = model.array(distance_warehouses)

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
route_quantity = model.sum(quantity[k][i] * model.contains(sequence, i) for i in range(nb_customers))
model.constraint(route_quantity <= truck_capacity)

# Distance traveled by each truck
dist_selector = model.lambda_function(lambda i:
model.at(distance_array, sequence[i - 1], sequence[i]))
route_distances[k] = model.sum(model.range(1, c), dist_selector) + \
model.iif(trucks_used[k], distance_warehouse_array[sequence] +
distance_warehouse_array[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[i])

total_distance = model.sum(route_distances)

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

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:
# 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]. +1 is to put it
# back in [1..nbCustomers]
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 warehouse 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_warehouses = compute_distance_warehouses(depot_x, depot_y, customers_x, customers_y)
return nb_customers, capacity, distance_matrix, distance_warehouses, demands

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

# Computes the distances to warehouse
def compute_distance_warehouses(depot_x, depot_y, customers_x, customers_y):
nb_customers = len(customers_x)
distance_warehouses = [None] * nb_customers
for i in range(nb_customers):
dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
distance_warehouses[i] = dist
return distance_warehouses

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\localsolver110.lib
sdvrp instances\S51D4.sd
Compilation / Execution (Linux)
g++ sdvrp.cpp -I/opt/localsolver_11_0/include -llocalsolver110 -lpthread -o sdvrp
./sdvrp instances/S51D4.sd
```/********** sdvrp.cpp **********/

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

using namespace localsolver;
using namespace std;

class Sdvrp {
public:
// Solver
LocalSolver localsolver;

// Number of customers
int nbCustomers;

// Capacity of the trucks
int truckCapacity;

// Demand of each customer
vector<lsint> demands;

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

// Distance between customers and warehouse
vector<lsint> distanceWarehouses;

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

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

void solve(int limit) {
// Declares 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, demands[i]);
}
}

// All customers must be visited by the trucks

// Creates distance as an array to be able to access it with an "at" operator
LSExpression distanceArray = model.array();
for (int n = 0; n < nbCustomers; n++) {
distanceMatrix[n].end()));
}
LSExpression distanceWarehousesArray = model.array(distanceWarehouses.begin(),
distanceWarehouses.end());

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

for (int k = 0; k < nbTrucks; k++) {
LSExpression sequence = customersSequences[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 routeQuantity = model.sum();
for (int i = 0; i < nbCustomers; i++)
model.contains(sequence, i));
model.constraint(routeQuantity <= truckCapacity);

// Distance traveled by truck k
LSExpression distSelector = model.createLambdaFunction([&](LSExpression i)
{ return model.at(distanceArray, model.at(sequence, model.sub(i, 1)),
model.at(sequence, i)); });
routeDistances[k] = model.sum(model.sum(model.range(1, c), distSelector),
model.iif(trucksUsed[k], model.sum (model.at(distanceWarehousesArray,
model.at(sequence,0)),model.at(distanceWarehousesArray,
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 >= demands[i]);
}

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

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

model.close();

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

localsolver.solve();
}

// Writes the solution in a file with the following format:
//  - each line k contain 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]. +1 is to put it back in
// [1..nbCustomers]
LSCollection customersCollection = customersSequences[k].getCollectionValue();
for (lsint i = 0; i < customersCollection.count(); i++) {
outfile << customersCollection[i] + 1 << " ";
}
outfile << endl;
}
}

private:

void readInputSdvrp(const string& fileName) {
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);
demands.resize(nbCustomers);
int demand;
for (int n = 0; n<nbCustomers;n++) {
stream >> demand;
demands[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();
}

// Computes the distance matrix
void computeDistanceMatrix(int depotX, int depotY, const vector<int>& customersX, const vector<int>& customersY) {
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(customersX[i], customersX[j], customersY[i], customersY[j]);
distanceMatrix[i][j] = distance;
distanceMatrix[j][i] = distance;
}
}

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

static lsint 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
{
// Solver
LocalSolver localsolver;

// Number of customers
int nbCustomers;

// Capacity of the trucks
int truckCapacity;

// Demand of each customer
long[] demands;

// Distance matrix
long[][] distanceMatrix;

// Distances between customers and warehouse
long[] distanceWarehouses;

// Number of trucks
int nbTrucks;

// Decision variables

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

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

// Distance traveled by all the trucks
LSExpression totalDistance;

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

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

// Reads instance data
{
}

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

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

trucksUsed = new LSExpression[nbTrucks];
customersSequences = new LSExpression[nbTrucks];
routeDistances = 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, demands[i]);
}

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

// Creates distance as an array to be able to access it with an "at" operator
LSExpression distanceArray = model.Array(distanceMatrix);
LSExpression distanceWarehouseArray = model.Array(distanceWarehouses);

for (int k = 0; k < nbTrucks; k++)
{
LSExpression sequence = customersSequences[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 routeQuantity = model.Sum();
for (int i = 0; i < nbCustomers; i++)
model.Contains(sequence, i)));
model.Constraint(routeQuantity <= truckCapacity);

// Distance traveled by truck k
LSExpression distSelector = model.LambdaFunction(i =>
distanceArray[sequence[i - 1], sequence[i]]);
routeDistances[k] = model.Sum(model.Range(1, c), distSelector)
+ model.If(trucksUsed[k], distanceWarehouseArray[sequence]
+ distanceWarehouseArray[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 >= demands[i]);
}

totalDistance = model.Sum(routeDistances);

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

model.Close();

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

localsolver.Solve();
}

// Writes the solution in a file with the following format:
//  - Each line k contain 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]. +1 is to put it
// back in [1..nbCustomers]
LSCollection customersCollection = customersSequences[k].GetCollectionValue();
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);
}
}

private static string[] SplitInput(StreamReader input)
{
string line = input.ReadLine();
if (line == null) return new string;
return line.Split(new[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
}

private void ReadInputSdvrp(string fileName)
{
{
string[] splitted;
splitted = SplitInput(input);

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

demands = new long[nbCustomers];

splitted = SplitInput(input);
for (int n = 0; n < nbCustomers; n++)
demands[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);
}
}

// Computes the distance matrix.
private void ComputeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY)
{
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(customersX[i], customersX[j], customersY[i], customersY[j]);
distanceMatrix[i][j] = dist;
distanceMatrix[j][i] = dist;
}
}
distanceWarehouses = new long[nbCustomers];
for (int i = 0; i < nbCustomers; ++i)
distanceWarehouses[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_11_0/bin/localsolver.jar
java -cp /opt/localsolver_11_0/bin/localsolver.jar:. Sdvrp instances/S51D4.sd
```/********** Sdvrp.java **********/
import java.util.*;
import java.io.*;
import java.lang.*;
import localsolver.*;

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

// Number of customers
int nbCustomers;

// Capacity of the trucks
private int truckCapacity;

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

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

// Distances between customers and warehouse
private long[] distanceWarehouses;

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

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

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

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

trucksUsed = new LSExpression[nbTrucks];
customersSequences = new LSExpression[nbTrucks];
routeDistances = 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, demands[i]);

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

// Creates distance as an array to be able to acces it with an "at" operator
LSExpression distanceArray = model.array(distanceMatrix);
LSExpression distanceWarehouseArray = model.array(distanceWarehouses);

for (int k = 0; k < nbTrucks; k++) {
LSExpression sequence = customersSequences[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 routeQuantity = model.sum();
for (int i = 0; i < nbCustomers; i++)
model.contains(sequence, i)));
model.constraint(model.leq(routeQuantity, truckCapacity));

// Distance traveled by truck k
LSExpression distSelector = model.lambdaFunction(i -> model.at(
distanceArray,
model.at(sequence, model.sub(i, 1)),
model.at(sequence, i)));
routeDistances[k] = model.sum(model.sum(model.range(1, c), distSelector),
model.iif(trucksUsed[k], model.sum(
model.at(distanceWarehouseArray, model.at(sequence, 0)),
model.at(distanceWarehouseArray, 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, demands[i]));
}

totalDistance = model.sum(routeDistances);

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

model.close();

// Parameterizes 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();
demands = new long[nbCustomers];
for (int n=1; n <= nbCustomers; n++) {
demands[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);
}
}

// Computes the distance matrix
private void computeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY) {
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(customersX[i], customersX[j], customersY[i], customersY[j]);
distanceMatrix[i][j] = dist;
distanceMatrix[j][i] = dist;
}
}

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

// Writes the solution in a file with the following format:
//  - Each line k contain 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]. +1 is to put it
// back in [1..nbCustomers]
LSCollection customersCollection = customersSequences[k].getCollectionValue();
for (int i = 0; i < customersCollection.count(); i++)
output.print((customersCollection.get(i) + 1) + " ");
output.println();
}
}
}
}
```