LocalSolver logo
is now
Hexaly logo

We're excited to share that we are moving forward. We're leaving behind the LocalSolver brand and transitioning to our new identity: Hexaly. This represents a leap forward in our mission to enable every organization to make better decisions faster when faced with operational and strategic challenges.

  • Docs »
  • Example tour »
  • Time Dependent Capacitated Vehicle Routing Problem with Time Windows (TDCVRPTW)

Time Dependent Capacitated Vehicle Routing Problem with Time Windows (TDCVRPTW)¶

Principles learned¶

  • Add multiple list decision variables

  • Use lambda expressions to define a recursive array

  • Define a sequence of expressions

Problem¶

../_images/tdcvrptw.svg

In the time dependent capacitated vehicle routing problem with time-windows (TDCVRPTW), a fleet of delivery vehicles with uniform capacity must service customers with known demand and opening hours for a single commodity. The vehicles start and end their routes at a common depot. Each customer can only be served by one vehicle. To model real-world traffic information, the travel time between customers depends on when the travel happens. The objectives are to minimize the fleet size and assign a sequence of customers to each truck of the fleet minimizing the total distance traveled such that all customers are served and the total demand served by each truck does not exceed its capacity.

In this example, the time-dependency is modeled using multiple constant travel time matrices. The horizon is divided into several time ranges. Each time range has its own travel time matrix for travel starting during this range.

Download the example


Data¶

The instances provided are based on the Solomon CVRPTW instances.

The format of the data files is as follows:

  • The first line gives the name of the instance

  • The fifth line contains the number of vehicles and their common capacity

  • From the 10th line, each line contains the integer data associated to each customer, starting with the depot:

    • The index of the customer

    • The x coordinate

    • The y coordinate

    • The demand

    • The earliest arrival

    • The latest arrival

    • The service time

To create TDCVRPTW instances from these CVRPTW instances, we chose to discretize the time horizon into five parts: early morning, morning peak, day, evening peak, and night. For each day part, the travel time between customers is modified according to coefficients depending on distance profiles. The shorter the distance is, the more the travel time is increased. This is a simple way to create non-proportionnal travel time matrices as it would be the case with real data.

Program¶

The Hexaly model is an extension of the CVRPTW model. We refer the reader to this model for the routing with time windows aspects of the problem.

To model the time-dependency, we introduce a timeToMatrixIdx array which gives for every time step of the horizon the index of the matrix that should be used.

We also have to define a travel time array with 3 dimensions. The first 2 dimensions represent the origin and destination customers. The third one is the index of the matrix. Hence, travelTime[origin][destination][matrixIdx] gives the travel time from origin to destination for the matrix at index matrixIdx.

Following the same logic, we define an array that gives the travel time to reach the depot depending on the matrix index. travelTimeWarehouse[customer][matrixIdx] is the travel time between customer and the depot for the matrix at index matrixIdx.

We can then modify the computation of the end times and home lateness defined in the classical CVRPTW. At each step and after doing the last delivery, we use the travel time given by the matrix valid at this time.

Execution:
localsolver tdcvrptw.lsp inFileName=instances/R101.25.txt [lsTimeLimit=] [solFileName=]
use io;

/* Read instance data. The input files follow the "Solomon" format*/
function input() {
    usage = "Usage: localsolver tdcvrptw.lsp "
            + "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;

    readInputCvrptw();

    computeDistanceMatrices();
}

/* Declare the optimization model */
function model() {
    customersSequences[k in 0...nbTrucks] <- list(nbCustomers);

    // All customers must be visited by exactly one truck
    constraint partition[k in 0...nbTrucks](customersSequences[k]);

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

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

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

        // End of each visit according to the traffic
        endTime[k] <- array(0...c, (i, prev) => max(earliestStart[sequence[i]],
            i == 0 ? travelTimeWarehouse[sequence[0]][timeToMatrixIdx[0]] :
            prev + travelTime[sequence[i-1]][sequence[i]][timeToMatrixIdx[round(prev)]])
            + serviceTime[sequence[i]]);

        // Arriving home after max horizon
        homeLateness[k] <- truckUsed[k] ? max(0, endTime[k][c-1] +
            travelTimeWarehouse[sequence[c-1]][timeToMatrixIdx[round(endTime[k][c-1])]] -
            maxHorizon) : 0;

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

        // Completing visit after latest end
        lateness[k] <- homeLateness[k] + sum(0...c,
                i => max(0, endTime[k][i] - latestEnd[sequence[i]]));
    }

    // Total lateness, must be 0 for a solution to be valid
    totalLateness <- sum[k in 0...nbTrucks](lateness[k]);

    // Total number of trucks used
    nbTrucksUsed <- sum[k in 0...nbTrucks](truckUsed[k]);

    // Total distance traveled (convention in Solomon's instances is to round to 2 decimals)
    totalDistance <- round(100 * sum[k in 0...nbTrucks](routeDistances[k])) / 100;

    // Objective: minimize the lateness, then the number of trucks used, then the distance traveled
    minimize totalLateness;
    minimize nbTrucksUsed;
    minimize totalDistance;
}


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

/* Write the solution in a file with the following format:
 * - number of trucks used and total distance
 * - for each truck the customers visited (omitting the start/end at the depot) */
function output() {
    if (solFileName == nil) return;
    local outfile = io.openWrite(solFileName);

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

function readInputCvrptw() {
    local inFile = io.openRead(inFileName);
    skipLines(inFile, 4);

    // Truck related data
    nbTrucks = inFile.readInt();
    truckCapacity = inFile.readInt();

    skipLines(inFile, 3);

    // Depot data
    local line = inFile.readln().split();
    depotIndex = line[0].toInt();
    depotX = line[1].toInt();
    depotY = line[2].toInt();
    maxHorizon = line[5].toInt();

    // Customers data
    i = 0;
    while (!inFile.eof()) {
        inLine = inFile.readln();
        line = inLine.split();
        if (count(line) == 0) break;
        if (count(line) != 7) throw "Wrong file format";
        customerIndex[i] = line[0].toInt();
        customerX[i] = line[1].toInt();
        customerY[i] = line[2].toInt();
        demands[i] = line[3].toInt();
        serviceTime[i] = line[6].toInt();
        earliestStart[i] = line[4].toInt();
        // in input files due date is meant as latest start time
        latestEnd[i] = line[5].toInt() + serviceTime[i];
        i = i + 1;
    }
    nbCustomers = i;

    inFile.close();
}

function skipLines(inFile, nbLines) {
    for [i in 0...nbLines] inFile.readln();
}

// Compute the distance matrices
function computeDistanceMatrices() {
    shortDistanceTravelTimeProfile = {1.00, 2.50, 1.75, 2.50, 1.00};
    mediumDistanceTravelTimeProfile = {1.00, 2.00, 1.50, 2.00, 1.00};
    longDistanceTravelTimeProfile = {1.00, 1.60, 1.10, 1.60, 1.00};
    travelTimeProfileMatrix = {
        shortDistanceTravelTimeProfile,
        mediumDistanceTravelTimeProfile,
        longDistanceTravelTimeProfile
    };
    distanceLevels = {10, 25};
    timeIntervalSteps = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
    nbTimeIntervals = count(timeIntervalSteps)-1;
    nbDistanceLevels = count(distanceLevels);

    for [i in 0...nbCustomers] {
        distanceMatrix[i][i] = 0;

        for [k in 0...nbTimeIntervals] {
            travelTime[i][i][k] = 0;
        }

        for [j in i+1...nbCustomers] {
            local localDistance = computeDist(i, j);
            distanceMatrix[j][i] = localDistance;
            distanceMatrix[i][j] = localDistance;

            local profileIdx = getProfile(localDistance);
            for [k in 0...nbTimeIntervals] {
                local localTravelTime = travelTimeProfileMatrix[profileIdx][k] * localDistance;
                travelTime[i][j][k] = localTravelTime;
                travelTime[j][i][k] = localTravelTime;
            }
        }
    }

    for [i in 0...nbTimeIntervals] {
        local timeStepStart = round(timeIntervalSteps[i] * maxHorizon);
        local timeStepEnd = round(timeIntervalSteps[i+1] * maxHorizon);
        for [j in timeStepStart...timeStepEnd+1] {
            timeToMatrixIdx[j] = i;
        }
    }

    for [i in 0...nbCustomers] {
        local localDistance = computeDepotDist(i);
        distanceDepot[i] = localDistance;

        local profileIdx = getProfile(localDistance);
        for [j in 0...nbTimeIntervals] {
            local localTravelTimeWareHouse = travelTimeProfileMatrix[profileIdx][j] * localDistance;
            travelTimeWarehouse[i][j] = localTravelTimeWareHouse;
        }
    }
}

function computeDist(i, j) {
    local x1 = customerX[i];
    local x2 = customerX[j];
    local y1 = customerY[i];
    local y2 = customerY[j];
    return computeDistance(x1, x2, y1, y2);
}

function computeDepotDist(i) {
    local x1 = customerX[i];
    local xd = depotX;
    local y1 = customerY[i];
    local yd = depotY;
    return computeDistance(x1, xd, y1, yd);
}

function computeDistance(x1, x2, y1, y2) {
    return sqrt(pow((x1 - x2), 2) + pow((y1 - y2), 2));
}

function getProfile(distance) {
    local idx = 0;
    while (idx < nbDistanceLevels && distance > distanceLevels[idx]) {
        idx = idx + 1;
    }
    return idx;
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python tdcvrptw.py instances\R101.25.txt
Execution (Linux)
export PYTHONPATH=/opt/localsolver_12_5/bin/python
python tdcvrptw.py instances/R101.25.txt
import localsolver
import sys
import math


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


def main(instance_file, str_time_limit, output_file):
    #
    # Read instance data
    #
    nb_customers, nb_trucks, truck_capacity, dist_matrix_data, travel_time_data, \
        time_to_matrix_idx_data, dist_depot_data, travel_time_warehouse_data,\
        demands_data, service_time_data, earliest_start_data, latest_end_data, \
        max_horizon = read_input_cvrptw(instance_file)

    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 k in range(nb_trucks)]

        # All customers must be visited by exactly one truck
        model.constraint(model.partition(customers_sequences))

        # Create LocalSolver arrays to be able to access them with an "at" operator
        demands = model.array(demands_data)
        earliest = model.array(earliest_start_data)
        latest = model.array(latest_end_data)
        service_time = model.array(service_time_data)
        dist_matrix = model.array()
        for n in range(nb_customers):
            dist_matrix.add_operand(model.array(dist_matrix_data[n]))
        travel_time = model.array()
        for n in range(nb_customers):
            time_matrix = model.array()
            for m in range(nb_customers):
                time_matrix.add_operand(model.array(travel_time_data[n][m]))
            travel_time.add_operand(time_matrix)
        time_to_matrix_idx = model.array(time_to_matrix_idx_data)
        dist_depot = model.array(dist_depot_data)
        travel_time_warehouse = model.array()
        for n in range(nb_customers):
            travel_time_warehouse.add_operand(model.array(travel_time_warehouse_data[n]))

        dist_routes = [None] * nb_trucks
        end_time = [None] * nb_trucks
        home_lateness = [None] * nb_trucks
        lateness = [None] * nb_trucks

        # A truck is used if it visits at least one customer
        trucks_used = [(model.count(customers_sequences[k]) > 0) for k in range(nb_trucks)]
        nb_trucks_used = model.sum(trucks_used)

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

            # The quantity needed in each route must not exceed the truck capacity
            demand_lambda = model.lambda_function(lambda j: demands[j])
            route_quantity = model.sum(sequence, demand_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]))
            dist_routes[k] = model.sum(model.range(1, c), dist_lambda) \
                + model.iif(c > 0, dist_depot[sequence[0]] + dist_depot[sequence[c - 1]], 0)

            # End of each visit according to the traffic
            end_time_lambda = model.lambda_function(
                lambda i, prev:
                    model.max(
                        earliest[sequence[i]],
                        model.iif(
                            i == 0, model.at(travel_time_warehouse, sequence[0], time_to_matrix_idx[0]),
                            prev + model.at(travel_time, sequence[i - 1], sequence[i],
                                            time_to_matrix_idx[model.round(prev)])))
                    + service_time[sequence[i]])

            end_time[k] = model.array(model.range(0, c), end_time_lambda)

            # Arriving home after max horizon
            home_lateness[k] = model.iif(
                trucks_used[k],
                model.max(
                    0,
                    end_time[k][c - 1] + model.at(travel_time_warehouse, sequence[c - 1],
                                                  time_to_matrix_idx[model.round(end_time[k][c - 1])]) - max_horizon),
                0)

            # Completing visit after latest end
            late_lambda = model.lambda_function(
                lambda i: model.max(0, end_time[k][i] - latest[sequence[i]]))
            lateness[k] = home_lateness[k] + model.sum(model.range(0, c), late_lambda)

        # Total lateness
        total_lateness = model.sum(lateness)

        # Total distance traveled
        total_distance = model.div(model.round(100 * model.sum(dist_routes)), 100)

        # Objective: minimize the number of trucks used, then minimize the distance traveled
        model.minimize(total_lateness)
        model.minimize(nb_trucks_used)
        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:
        #  - number of trucks used and total distance
        #  - for each truck the customers visited (omitting the start/end at the depot)
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("%d %d\n" % (nb_trucks_used.value, 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 as in the data files (0 being the depot)
                    for customer in customers_sequences[k].value:
                        f.write("%d " % (customer + 1))
                    f.write("\n")


# The input files follow the "Solomon" format
def read_input_cvrptw(filename):
    file_it = iter(read_elem(filename))

    for i in range(4):
        next(file_it)

    nb_trucks = int(next(file_it))
    truck_capacity = int(next(file_it))

    for i in range(13):
        next(file_it)

    depot_x = int(next(file_it))
    depot_y = int(next(file_it))

    for i in range(2):
        next(file_it)

    max_horizon = int(next(file_it))

    next(file_it)

    customers_x = []
    customers_y = []
    demands = []
    earliest_start = []
    latest_end = []
    service_time = []

    while True:
        val = next(file_it, None)
        if val is None:
            break
        i = int(val) - 1
        customers_x.append(int(next(file_it)))
        customers_y.append(int(next(file_it)))
        demands.append(int(next(file_it)))
        ready = int(next(file_it))
        due = int(next(file_it))
        stime = int(next(file_it))
        earliest_start.append(ready)
        # in input files due date is meant as latest start time
        latest_end.append(due + stime)
        service_time.append(stime)

    nb_customers = i + 1

    short_distance_travel_time_profile = [1.00, 2.50, 1.75, 2.50, 1.00]
    medium_distance_travel_time_profile = [1.00, 2.00, 1.50, 2.00, 1.00]
    long_distance_travel_time_profile = [1.00, 1.60, 1.10, 1.60, 1.00]
    travel_time_profile_matrix = [
        short_distance_travel_time_profile,
        medium_distance_travel_time_profile,
        long_distance_travel_time_profile
    ]
    distance_levels = [10, 25]
    time_interval_steps = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
    nb_time_intervals = len(time_interval_steps) - 1
    nb_distance_levels = len(distance_levels)

    # Compute distance matrices
    distance_matrix, travel_time, time_to_matrix_idx = compute_distance_matrices(
        customers_x, customers_y, max_horizon, travel_time_profile_matrix, time_interval_steps, nb_time_intervals,
        distance_levels, nb_distance_levels)
    distance_depots, travel_time_warehouse = compute_distance_depots(
        depot_x, depot_y, customers_x, customers_y, travel_time_profile_matrix, nb_time_intervals, distance_levels,
        nb_distance_levels)

    return nb_customers, nb_trucks, truck_capacity, distance_matrix, travel_time, time_to_matrix_idx, distance_depots, \
        travel_time_warehouse, demands, service_time, earliest_start, latest_end, max_horizon


# Computes the distance matrices
def compute_distance_matrices(customers_x, customers_y, max_horizon, travel_time_profile_matrix,
                              time_interval_steps, nb_time_intervals, distance_levels, nb_distance_levels):

    nb_customers = len(customers_x)
    distance_matrix = [[None for _ in range(nb_customers)] for _ in range(nb_customers)]
    travel_time = [[[None for _ in range(nb_time_intervals)] for _ in range(nb_customers)] for _ in range(nb_customers)]
    time_to_matrix_idx = [None for _ in range(max_horizon)]
    for i in range(nb_customers):
        distance_matrix[i][i] = 0

        for k in range(nb_time_intervals):
            travel_time[i][i][k] = 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

            profile_idx = get_profile(dist, distance_levels, nb_distance_levels)
            for k in range(nb_time_intervals):
                local_travel_time = travel_time_profile_matrix[profile_idx][k] * dist
                travel_time[i][j][k] = local_travel_time
                travel_time[j][i][k] = local_travel_time

    for i in range(nb_time_intervals):
        time_step_start = int(round(time_interval_steps[i] * max_horizon))
        time_step_end = int(round(time_interval_steps[i+1] * max_horizon))
        for j in range(time_step_start, time_step_end):
            time_to_matrix_idx[j] = i
    return distance_matrix, travel_time, time_to_matrix_idx


# Computes the distances to depot
def compute_distance_depots(depot_x, depot_y, customers_x, customers_y, travel_time_profile_matrix,
                            nb_time_intervals, distance_levels, nb_distance_levels):
    nb_customers = len(customers_x)
    distance_depots = [None] * nb_customers
    travel_time_warehouse = [[None for _ in range(nb_time_intervals)] for _ in range(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

        profile_idx = get_profile(dist, distance_levels, nb_distance_levels)
        for j in range(nb_time_intervals):
            local_travel_time_warehouse = travel_time_profile_matrix[profile_idx][j] * dist
            travel_time_warehouse[i][j] = local_travel_time_warehouse
    return distance_depots, travel_time_warehouse


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


def get_profile(dist, distance_levels, nb_distance_levels):
    idx = 0
    while idx < nb_distance_levels and dist > distance_levels[idx]:
        idx += 1
    return idx


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

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

    main(instance_file, str_time_limit, output_file)
Compilation / Execution (Windows)
cl /EHsc tdcvrptw.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver125.lib
tdcvrptw instances\R101.25.txt
Compilation / Execution (Linux)
g++ tdcvrptw.cpp -I/opt/localsolver_12_5/include -llocalsolver125 -lpthread -o tdcvrptw
./tdcvrptw instances/R101.25.txt
#include "localsolver.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace localsolver;
using namespace std;

class Tdcvrptw {
public:
    // LocalSolver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Latest allowed arrival to depot
    int maxHorizon;

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

    // Earliest arrival for each customer
    vector<int> earliestStartData;

    // Latest departure from each customer
    vector<int> latestEndData;

    // Service time for each customer
    vector<int> serviceTimeData;

    // Distance matrix between customers
    vector<vector<double>> distMatrixData;

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

    // Travel time coefficients for each profile
    vector<double> shortDistanceTravelTimeProfile{1.00, 2.50, 1.75, 2.50, 1.00};
    vector<double> mediumDistanceTravelTimeProfile{1.00, 2.00, 1.50, 2.00, 1.00};
    vector<double> longDistanceTravelTimeProfile{1.00, 1.60, 1.10, 1.60, 1.00};
    vector<vector<double>> travelTimeProfileMatrix{shortDistanceTravelTimeProfile, mediumDistanceTravelTimeProfile,
                                                   longDistanceTravelTimeProfile};

    // Distance levels
    vector<int> distanceLevels{10, 25};

    // Intervals of the temporal discretization
    vector<double> timeIntervalSteps{0.0, 0.2, 0.4, 0.6, 0.8, 1.0};

    // Number of time intervals
    int nbTimeIntervals = timeIntervalSteps.size() - 1;

    // Number of distance levels
    int nbDistanceLevels = distanceLevels.size();

    // Travel time between customers for each day part
    vector<vector<vector<double>>> travelTimeData;

    // Time interval index for each time unit
    vector<int> timeToMatrixIdxData;

    // Travel time between customers and depot for each day part
    vector<vector<double>> travelTimeWarehouseData;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    vector<LSExpression> customersSequences;

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

    // Cumulated lateness in the solution (must be 0 for the solution to be valid)
    LSExpression totalLateness;

    // Number of trucks used in the solution
    LSExpression nbTrucksUsed;

    // Distance traveled by all the trucks
    LSExpression totalDistance;

    Tdcvrptw() {}

    /* Read instance data */
    void readInstance(const string& fileName) { readInputCvrptw(fileName); }

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

        // Sequence of customers visited by each truck
        customersSequences.resize(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            customersSequences[k] = model.listVar(nbCustomers);
        }

        // All customers must be visited by exactly one truck
        model.constraint(model.partition(customersSequences.begin(), customersSequences.end()));

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = model.array(demandsData.begin(), demandsData.end());
        LSExpression earliest = model.array(earliestStartData.begin(), earliestStartData.end());
        LSExpression latest = model.array(latestEndData.begin(), latestEndData.end());
        LSExpression serviceTime = model.array(serviceTimeData.begin(), serviceTimeData.end());
        LSExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            distMatrix.addOperand(model.array(distMatrixData[n].begin(), distMatrixData[n].end()));
        }
        LSExpression travelTime = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            LSExpression timeMatrix = model.array();
            for (int m = 0; m < nbCustomers; ++m) {
                timeMatrix.addOperand(model.array(travelTimeData[n][m].begin(), travelTimeData[n][m].end()));
            }
            travelTime.addOperand(timeMatrix);
        }
        LSExpression timeToMatrixIdx = model.array(timeToMatrixIdxData.begin(), timeToMatrixIdxData.end());
        LSExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());
        LSExpression travelTimeWarehouse = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            travelTimeWarehouse.addOperand(
                model.array(travelTimeWarehouseData[n].begin(), travelTimeWarehouseData[n].end()));
        }

        trucksUsed.resize(nbTrucks);
        vector<LSExpression> distRoutes(nbTrucks), endTime(nbTrucks), homeLateness(nbTrucks), lateness(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 needed in each route must not exceed the truck capacity
            LSExpression demandLambda =
                model.createLambdaFunction([&](LSExpression j) { return demands[j]; });
            LSExpression routeQuantity = model.sum(sequence, demandLambda);
            model.constraint(routeQuantity <= truckCapacity);

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

            // End of each visit according to the traffic
            LSExpression endTimeLambda = model.createLambdaFunction([&](LSExpression i, LSExpression prev) {
                return model.max(earliest[sequence[i]],
                                 model.iif(i == 0, model.at(travelTimeWarehouse, sequence[0], timeToMatrixIdx[0]),
                                           prev + model.at(travelTime, sequence[i - 1], sequence[i],
                                                           timeToMatrixIdx[model.round(prev)]))) +
                       serviceTime[sequence[i]];
            });

            endTime[k] = model.array(model.range(0, c), endTimeLambda);

            // Arriving home after max horizon
            homeLateness[k] = model.iif(trucksUsed[k],
                                        model.max(0, endTime[k][c - 1] +
                                                         model.at(travelTimeWarehouse, sequence[c - 1],
                                                                  timeToMatrixIdx[model.round(endTime[k][c - 1])]) -
                                                         maxHorizon),
                                        0);

            // Completing visit after latest end
            LSExpression lateLambda = model.createLambdaFunction(
                [&](LSExpression i) { return model.max(0, endTime[k][i] - latest[sequence[i]]); });
            lateness[k] = homeLateness[k] + model.sum(model.range(0, c), lateLambda);
        }

        // Total lateness
        totalLateness = model.sum(lateness.begin(), lateness.end());

        // Total number of trucks used
        nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());

        // Total distance traveled (convention in Solomon's instances is to round to 2 decimals)
        totalDistance = model.round(100 * model.sum(distRoutes.begin(), distRoutes.end())) / 100;

        // Objective: minimize the lateness, then the number of trucks used, then the distance traveled
        model.minimize(totalLateness);
        model.minimize(nbTrucksUsed);
        model.minimize(totalDistance);

        model.close();

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

        localsolver.solve();
    }

    /* Write the solution in a file with the following format:
     *  - number of trucks used and total distance
     *  - for each truck the customers visited (omitting the start/end at the depot) */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << nbTrucksUsed.getValue() << " " << totalDistance.getDoubleValue() << 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
            // as in the data files (0 being the depot)
            LSCollection customersCollection = customersSequences[k].getCollectionValue();
            for (int i = 0; i < customersCollection.count(); ++i) {
                outfile << customersCollection[i] + 1 << " ";
            }
            outfile << endl;
        }
    }

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

        string str;
        long tmp;

        int depotX, depotY;
        vector<int> customersX;
        vector<int> customersY;

        getline(infile, str);
        getline(infile, str);
        getline(infile, str);
        getline(infile, str);

        infile >> nbTrucks;
        infile >> truckCapacity;

        getline(infile, str);
        getline(infile, str);
        getline(infile, str);
        getline(infile, str);

        infile >> tmp;
        infile >> depotX;
        infile >> depotY;
        infile >> tmp;
        infile >> tmp;
        infile >> maxHorizon;
        infile >> tmp;

        while (infile >> tmp) {
            int cx, cy, demand, ready, due, service;
            infile >> cx;
            infile >> cy;
            infile >> demand;
            infile >> ready;
            infile >> due;
            infile >> service;

            customersX.push_back(cx);
            customersY.push_back(cy);
            demandsData.push_back(demand);
            earliestStartData.push_back(ready);
            latestEndData.push_back(due + service); // in input files due date is meant as latest start time
            serviceTimeData.push_back(service);
        }

        nbCustomers = customersX.size();

        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);
        travelTimeData.resize(nbCustomers);
        timeToMatrixIdxData.resize(maxHorizon);
        travelTimeWarehouseData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i].resize(nbCustomers);
            travelTimeData[i].resize(nbCustomers);
            for (int j = 0; j < nbCustomers; ++j) {
                travelTimeData[i][j].resize(nbTimeIntervals);
            }
        }
        for (int i = 0; i < nbCustomers; ++i) {
            travelTimeWarehouseData[i].resize(nbTimeIntervals);
        }
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int k = 0; k < nbTimeIntervals; ++k) {
                travelTimeData[i][i][k] = 0;
            }
            for (int j = i + 1; j < nbCustomers; ++j) {
                double distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = distance;
                distMatrixData[j][i] = distance;

                int profileIdx = getProfile(distance);
                for (int k = 0; k < nbTimeIntervals; ++k) {
                    double localTravelTime = travelTimeProfileMatrix[profileIdx][k] * distance;
                    travelTimeData[i][j][k] = localTravelTime;
                    travelTimeData[j][i][k] = localTravelTime;
                }
            }
        }

        for (int i = 0; i < nbTimeIntervals; ++i) {
            int timeStepStart = (int)round(timeIntervalSteps[i] * maxHorizon);
            int timeStepEnd = (int)round(timeIntervalSteps[i + 1] * maxHorizon);
            for (int j = timeStepStart; j < timeStepEnd; ++j) {
                timeToMatrixIdxData[j] = i;
            }
        }

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

            int profileIdx = getProfile(distance);
            for (int j = 0; j < nbTimeIntervals; ++j) {
                double localTravelTimeWarehouse = travelTimeProfileMatrix[profileIdx][j] * distance;
                travelTimeWarehouseData[i][j] = localTravelTimeWarehouse;
            }
        }
    }

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

    int getProfile(double distance) {
        int idx = 0;
        while (idx < nbDistanceLevels && distance > distanceLevels[idx]) {
            idx += 1;
        }
        return idx;
    }
};

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

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

    try {
        Tdcvrptw model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc Tdcvrptw.cs /reference:localsolvernet.dll
Tdcvrptw instances\R101.25.txt
using System;
using System.IO;
using System.Collections.Generic;
using localsolver;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Latest allowed arrival to depot
    int maxHorizon;

    // Demand for each customer
    List<int> demandsData;

    // Earliest arrival for each customer
    List<int> earliestStartData;

    // Latest departure from each customer
    List<int> latestEndData;

    // Service time for each customer
    List<int> serviceTimeData;

    // Distance matrix between customers
    double[][] distMatrixData;

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

    // Travel time coefficients for each profile
    static double[] shortDistanceTravelTimeProfile = { 1.00, 2.50, 1.75, 2.50, 1.00 };
    static double[] mediumDistanceTravelTimeProfile = { 1.00, 2.00, 1.50, 2.00, 1.00 };
    static double[] longDistanceTravelTimeProfile = { 1.00, 1.60, 1.10, 1.60, 1.00 };
    static double[][] travelTimeProfileMatrix =
    {
        shortDistanceTravelTimeProfile,
        mediumDistanceTravelTimeProfile,
        longDistanceTravelTimeProfile
    };

    // Distance levels
    static int[] distanceLevels = { 10, 25 };

    // Intervals of the temporal discretization
    static double[] timeIntervalSteps = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };

    // Number of time intervals
    static int nbTimeIntervals = timeIntervalSteps.Length - 1;

    // Number of distances levels
    static int nbDistanceLevels = distanceLevels.Length;

    // Travel time between customers for each day part
    double[][][] travelTimeData;

    // Time interval index for each time unit
    int[] timeToMatrixIdxData;

    // Travel time between customers and depot for each day part
    double[][] travelTimeWarehouseData;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    LSExpression[] customersSequences;

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

    // Cumulated lateness in the solution (must be 0 for the solution to be valid)
    LSExpression totalLateness;

    // Number of trucks used in the solution
    LSExpression nbTrucksUsed;

    // Distance traveled by all the trucks
    LSExpression totalDistance;

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

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputCvrptw(fileName);
    }

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

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

        trucksUsed = new LSExpression[nbTrucks];
        customersSequences = new LSExpression[nbTrucks];
        LSExpression[] distRoutes = new LSExpression[nbTrucks];
        LSExpression[] endTime = new LSExpression[nbTrucks];
        LSExpression[] homeLateness = new LSExpression[nbTrucks];
        LSExpression[] lateness = new LSExpression[nbTrucks];

        // Sequence of customers visited by each truck
        for (int k = 0; k < nbTrucks; ++k)
            customersSequences[k] = model.List(nbCustomers);

        // All customers must be visited by exactly one truck
        model.Constraint(model.Partition(customersSequences));

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = model.Array(demandsData);
        LSExpression earliest = model.Array(earliestStartData);
        LSExpression latest = model.Array(latestEndData);
        LSExpression serviceTime = model.Array(serviceTimeData);
        LSExpression distDepot = model.Array(distDepotData);
        LSExpression distMatrix = model.Array(distMatrixData);
        LSExpression travelTime = model.Array(travelTimeData);
        LSExpression timeToMatrixIdx = model.Array(timeToMatrixIdxData);
        LSExpression travelTimeWarehouse = model.Array(travelTimeWarehouseData);

        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 needed in each route must not exceed the truck capacity
            LSExpression demandLambda = model.LambdaFunction(j => demands[j]);
            LSExpression routeQuantity = model.Sum(sequence, demandLambda);
            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(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

            // End of each visit according to the traffic
            LSExpression endTimeLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        earliest[sequence[i]],
                        model.If(
                            i == 0,
                            travelTimeWarehouse[sequence[0], timeToMatrixIdx[0]],
                            prev
                                + model.At(
                                    travelTime,
                                    sequence[i - 1],
                                    sequence[i],
                                    timeToMatrixIdx[model.Round(prev)]
                                )
                        )
                    ) + serviceTime[sequence[i]]
            );

            endTime[k] = model.Array(model.Range(0, c), endTimeLambda);

            // Arriving home after max_horizon
            homeLateness[k] = model.If(
                trucksUsed[k],
                model.Max(
                    0,
                    endTime[k][c - 1]
                        + travelTimeWarehouse[
                            sequence[c - 1],
                            timeToMatrixIdx[model.Round(endTime[k][c - 1])]
                        ]
                        - maxHorizon
                ),
                0
            );

            // Completing visit after latest_end
            LSExpression lateLambda = model.LambdaFunction(
                i => model.Max(endTime[k][i] - latest[sequence[i]], 0)
            );
            lateness[k] = homeLateness[k] + model.Sum(model.Range(0, c), lateLambda);
        }

        // Total lateness
        totalLateness = model.Sum(lateness);

        // Total number of trucks used
        nbTrucksUsed = model.Sum(trucksUsed);

        // Total distance traveled (convention in Solomon's instances is to round to 2 decimals)
        totalDistance = model.Round(100 * model.Sum(distRoutes)) / 100;

        // Objective: minimize the lateness, then the number of trucks used, then the distance traveled
        model.Minimize(totalLateness);
        model.Minimize(nbTrucksUsed);
        model.Minimize(totalDistance);

        model.Close();

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

        localsolver.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - number of trucks used and total distance
     *  - for each truck the customers visited (omitting the start/end at the depot) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(nbTrucksUsed.GetValue() + " " + totalDistance.GetDoubleValue());
            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
                // as in the data files (0 being the depot)
                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: Tdcvrptw inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

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

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

    // The input files follow the "Solomon" format
    void ReadInputCvrptw(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;

            input.ReadLine();
            input.ReadLine();
            input.ReadLine();
            input.ReadLine();

            splitted = SplitInput(input);
            nbTrucks = int.Parse(splitted[0]);
            truckCapacity = int.Parse(splitted[1]);

            input.ReadLine();
            input.ReadLine();
            input.ReadLine();
            input.ReadLine();

            splitted = SplitInput(input);
            int depotX = int.Parse(splitted[1]);
            int depotY = int.Parse(splitted[2]);
            maxHorizon = int.Parse(splitted[5]);

            List<int> customersX = new List<int>();
            List<int> customersY = new List<int>();
            demandsData = new List<int>();
            earliestStartData = new List<int>();
            latestEndData = new List<int>();
            serviceTimeData = new List<int>();

            while (!input.EndOfStream)
            {
                splitted = SplitInput(input);
                if (splitted.Length < 7)
                    break;
                customersX.Add(int.Parse(splitted[1]));
                customersY.Add(int.Parse(splitted[2]));
                demandsData.Add(int.Parse(splitted[3]));
                int ready = int.Parse(splitted[4]);
                int due = int.Parse(splitted[5]);
                int service = int.Parse(splitted[6]);

                earliestStartData.Add(ready);
                latestEndData.Add(due + service); // in input files due date is meant as latest start time
                serviceTimeData.Add(service);
            }

            nbCustomers = customersX.Count;

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

    // Compute the distance matrix
    void ComputeDistanceMatrix(int depotX, int depotY, List<int> customersX, List<int> customersY)
    {
        distMatrixData = new double[nbCustomers][];
        travelTimeData = new double[nbCustomers][][];
        timeToMatrixIdxData = new int[maxHorizon];
        travelTimeWarehouseData = new double[nbCustomers][];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i] = new double[nbCustomers];
            double[][] timeMatrix = new double[nbCustomers][];
            for (int j = 0; j < nbCustomers; ++j)
                timeMatrix[j] = new double[nbTimeIntervals];
            travelTimeData[i] = timeMatrix;
        }

        for (int i = 0; i < nbCustomers; ++i)
            travelTimeWarehouseData[i] = new double[nbTimeIntervals];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int k = 0; k < nbTimeIntervals; ++k)
                travelTimeData[i][i][k] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                double dist = ComputeDist(
                    customersX[i],
                    customersX[j],
                    customersY[i],
                    customersY[j]
                );
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;

                int profileIdx = GetProfile(dist);
                for (int k = 0; k < nbTimeIntervals; ++k)
                {
                    double localTravelTime = travelTimeProfileMatrix[profileIdx][k] * dist;
                    travelTimeData[i][j][k] = localTravelTime;
                    travelTimeData[j][i][k] = localTravelTime;
                }
            }
        }

        for (int i = 0; i < nbTimeIntervals; ++i)
        {
            int timeStepStart = (int)Math.Round(timeIntervalSteps[i] * maxHorizon);
            int timeStepEnd = (int)Math.Round(timeIntervalSteps[i + 1] * maxHorizon);
            for (int j = timeStepStart; j < timeStepEnd; ++j)
                timeToMatrixIdxData[j] = i;
        }

        distDepotData = new double[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
        {
            double dist = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
            distDepotData[i] = dist;

            int profileIdx = GetProfile(dist);
            for (int j = 0; j < nbTimeIntervals; ++j)
            {
                double localTravelTimeWarehouse = travelTimeProfileMatrix[profileIdx][j] * dist;
                travelTimeWarehouseData[i][j] = localTravelTimeWarehouse;
            }
        }
    }

    double ComputeDist(int xi, int xj, int yi, int yj)
    {
        return Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
    }

    int GetProfile(double dist)
    {
        int idx = 0;
        while (idx < nbDistanceLevels && dist > distanceLevels[idx])
            idx += 1;
        return idx;
    }
}
Compilation / Execution (Windows)
javac Tdcvrptw.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Tdcvrptw instances\R101.25.txt
Compilation / Execution (Linux)
javac Tdcvrptw.java -cp /opt/localsolver_12_5/bin/localsolver.jar
java -cp /opt/localsolver_12_5/bin/localsolver.jar:. Tdcvrptw instances/R101.25.txt
import java.util.*;
import java.io.*;
import localsolver.*;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    private int truckCapacity;

    // Latest allowed arrival to depot
    int maxHorizon;

    // Demand for each customer
    List<Integer> demandsData;

    // Earliest arrival for each customer
    List<Integer> earliestStartData;

    // Latest departure from each customer
    List<Integer> latestEndData;

    // Service time for each customer
    List<Integer> serviceTimeData;

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

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

    // Travel time coefficients for each profile
    private static double shortDistanceTravelTimeProfile[] = { 1.00, 2.50, 1.75, 2.50, 1.00 };
    private static double mediumDistanceTravelTimeProfile[] = { 1.00, 2.00, 1.50, 2.00, 1.00 };
    private static double longDistanceTravelTimeProfile[] = { 1.00, 1.60, 1.10, 1.60, 1.00 };
    private static double[] travelTimeProfileMatrix[] = {
            shortDistanceTravelTimeProfile,
            mediumDistanceTravelTimeProfile,
            longDistanceTravelTimeProfile
    };

    // Distance levels
    private static int distanceLevels[] = { 10, 25 };

    // Intervals of the temporal discretization
    private static double timeIntervalSteps[] = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };

    // Number of time intervals
    private static int nbTimeIntervals = timeIntervalSteps.length - 1;

    // Number of distances levels
    private static int nbDistanceLevels = distanceLevels.length;

    // Travel time between customers for each day part
    private double[][][] travelTimeData;

    // Time interval index for each time unit
    private int[] timeToMatrixIdxData;

    // Travel time between customers and depot for each day part
    private double[][] travelTimeWarehouseData;

    // Number of trucks
    private int nbTrucks;

    // Decision variables
    private LSExpression[] customersSequences;

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

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

    // End time array for each truck
    private LSExpression[] endTime;

    // Home lateness for each truck
    private LSExpression[] homeLateness;

    // Cumulated Lateness for each truck
    private LSExpression[] lateness;

    // Cumulated lateness in the solution (must be 0 for the solution to be valid)
    private LSExpression totalLateness;

    // Number of trucks used in the solution
    private LSExpression nbTrucksUsed;

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

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

    /* Read instance data */
    private void readInstance(String fileName) throws IOException {
        readInputCvrptw(fileName);
    }

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

        trucksUsed = new LSExpression[nbTrucks];
        customersSequences = new LSExpression[nbTrucks];
        distRoutes = new LSExpression[nbTrucks];
        endTime = new LSExpression[nbTrucks];
        homeLateness = new LSExpression[nbTrucks];
        lateness = new LSExpression[nbTrucks];

        // Sequence of customers visited by each truck.
        for (int k = 0; k < nbTrucks; ++k)
            customersSequences[k] = m.listVar(nbCustomers);

        // All customers must be visited by exactly one truck
        m.constraint(m.partition(customersSequences));

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = m.array(demandsData);
        LSExpression earliest = m.array(earliestStartData);
        LSExpression latest = m.array(latestEndData);
        LSExpression serviceTime = m.array(serviceTimeData);
        LSExpression distDepot = m.array(distDepotData);
        LSExpression distMatrix = m.array(distMatrixData);
        LSExpression travelTime = m.array(travelTimeData);
        LSExpression timeToMatrixIdx = m.array(timeToMatrixIdxData);
        LSExpression travelTimeWarehouse = m.array(travelTimeWarehouseData);

        for (int k = 0; k < nbTrucks; ++k) {
            LSExpression sequence = customersSequences[k];
            LSExpression c = m.count(sequence);

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

            // The quantity needed in each route must not exceed the truck capacity
            LSExpression demandLambda = m.lambdaFunction(j -> m.at(demands, j));
            LSExpression routeQuantity = m.sum(sequence, demandLambda);
            m.constraint(m.leq(routeQuantity, truckCapacity));

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

            // End of each visit according to the traffic
            LSExpression endTimeLambda = m.lambdaFunction((i, prev) -> m.sum(
                    m.max(m.at(earliest, m.at(sequence, i)),
                            m.sum(m.iif(m.eq(i, 0),
                                    m.at(travelTimeWarehouse, m.at(sequence, 0), m.at(timeToMatrixIdx, 0)),
                                    m.sum(prev, m.at(travelTime, m.at(sequence, m.sub(i, 1)), m.at(sequence, i),
                                            m.at(timeToMatrixIdx, m.round(prev))))))),
                    m.at(serviceTime, m.at(sequence, i))));

            endTime[k] = m.array(m.range(0, c), endTimeLambda);

            LSExpression theEnd = endTime[k];

            // Arriving home after max_horizon
            homeLateness[k] = m.iif(trucksUsed[k],
                    m.max(0,
                            m.sum(m.at(theEnd, m.sub(c, 1)),
                                    m.sub(m.at(m.at(travelTimeWarehouse, m.at(sequence, m.sub(c, 1))),
                                            m.at(timeToMatrixIdx, m.round(m.at(theEnd, m.sub(c, 1))))), maxHorizon))),
                    0);

            // Completing visit after latest_end
            LSExpression lateLambda = m
                    .lambdaFunction(i -> m.max(m.sub(m.at(theEnd, i), m.at(latest, m.at(sequence, i))), 0));
            lateness[k] = m.sum(homeLateness[k], m.sum(m.range(0, c), lateLambda));
        }

        totalLateness = m.sum(lateness);
        nbTrucksUsed = m.sum(trucksUsed);
        totalDistance = m.div(m.round(m.prod(100, m.sum(distRoutes))), 100);

        // Objective: minimize the number of trucks used, then minimize the distance traveled
        m.minimize(totalLateness);
        m.minimize(nbTrucksUsed);
        m.minimize(totalDistance);

        m.close();

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

        localsolver.solve();
    }

    // Write the solution in a file with the following format:
    // - number of trucks used and total distance
    // - for each truck the customers visited (omitting the start/end at the depot)
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(nbTrucksUsed.getValue() + " " + totalDistance.getDoubleValue());
            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
                // as in the data files (0 being the depot)
                LSCollection customersCollection = customersSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); ++i) {
                    output.print((customersCollection.get(i) + 1) + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "Solomon" format
    private void readInputCvrptw(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.nextLine();
            input.nextLine();
            input.nextLine();
            input.nextLine();

            nbTrucks = input.nextInt();
            truckCapacity = input.nextInt();

            input.nextLine();
            input.nextLine();
            input.nextLine();
            input.nextLine();

            input.nextInt();
            int depotX = input.nextInt();
            int depotY = input.nextInt();
            input.nextInt();
            input.nextInt();
            maxHorizon = input.nextInt();
            input.nextInt();

            List<Integer> customersX = new ArrayList<Integer>();
            List<Integer> customersY = new ArrayList<Integer>();
            demandsData = new ArrayList<Integer>();
            earliestStartData = new ArrayList<Integer>();
            latestEndData = new ArrayList<Integer>();
            serviceTimeData = new ArrayList<Integer>();

            while (input.hasNextInt()) {
                input.nextInt();
                int cx = input.nextInt();
                int cy = input.nextInt();
                int demand = input.nextInt();
                int ready = input.nextInt();
                int due = input.nextInt();
                int service = input.nextInt();

                customersX.add(cx);
                customersY.add(cy);
                demandsData.add(demand);
                earliestStartData.add(ready);
                latestEndData.add(due + service);// in input files due date is meant as latest start time
                serviceTimeData.add(service);
            }

            nbCustomers = customersX.size();

            computeDistanceMatrix(depotX, depotY, customersX, customersY);

        }
    }

    // Computes the distance matrix
    private void computeDistanceMatrix(int depotX, int depotY, List<Integer> customersX, List<Integer> customersY) {
        distMatrixData = new double[nbCustomers][nbCustomers];
        travelTimeData = new double[nbCustomers][nbCustomers][nbTimeIntervals];
        timeToMatrixIdxData = new int[maxHorizon];
        travelTimeWarehouseData = new double[nbCustomers][nbTimeIntervals];
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int k = 0; k < nbTimeIntervals; ++k) {
                travelTimeData[i][i][k] = 0;
            }
            for (int j = i + 1; j < nbCustomers; ++j) {
                double dist = computeDist(customersX.get(i), customersX.get(j), customersY.get(i), customersY.get(j));
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;

                int profileIdx = getProfile(dist);
                for (int k = 0; k < nbTimeIntervals; ++k) {
                    double localTravelTime = travelTimeProfileMatrix[profileIdx][k] * dist;
                    travelTimeData[i][j][k] = localTravelTime;
                    travelTimeData[j][i][k] = localTravelTime;
                }
            }
        }

        for (int i = 0; i < nbTimeIntervals; ++i) {
            int timeStepStart = (int) Math.round(timeIntervalSteps[i] * maxHorizon);
            int timeStepEnd = (int) Math.round(timeIntervalSteps[i + 1] * maxHorizon);
            for (int j = timeStepStart; j < timeStepEnd; ++j) {
                timeToMatrixIdxData[j] = i;
            }
        }

        distDepotData = new double[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            double dist = computeDist(depotX, customersX.get(i), depotY, customersY.get(i));
            distDepotData[i] = dist;

            int profileIdx = getProfile(dist);
            for (int j = 0; j < nbTimeIntervals; ++j) {
                double localTravelTimeWarehouse = travelTimeProfileMatrix[profileIdx][j] * dist;
                travelTimeWarehouseData[i][j] = localTravelTimeWarehouse;
            }
        }
    }

    private double computeDist(int xi, int xj, int yi, int yj) {
        return Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
    }

    private int getProfile(double dist) {
        int idx = 0;
        while (idx < nbDistanceLevels && dist > distanceLevels[idx]) {
            idx += 1;
        }
        return idx;
    }

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

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

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