Pickup and Delivery with Time Windows (PDPTW)

Principles learned

  • Add multiple list decision variables

  • Use a lambda expression to define a recursive array

  • Define a sequence of expressions

  • Use of “contains” and “indexOf” operators

Problem

../_images/pdp.svg

In the pickup-and-delivery problem with time-windows (PDPTW), a fleet of delivery vehicles with uniform capacity must collect and deliver items according to the demand of customers, and must do so during their opening hours.

The vehicles start and end their routes at a common depot. Each customer can only be served by one vehicle. The objectives are to minimize the fleet size and to assign a sequence of customers to each truck of the fleet minimizing the total distance traveled such that each item is picked up then delivered by the same vehicle, and the capacity of the truck is respected at any point in the tour.

Download the example


Data

The instances provided come from the Li & Lim instances.

The format of the data files is as follows:

  • The first line gives the number of vehicles, the capacity and the speed (not used)

  • From the 2nd 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

    • The index of the corresponding pickup order (0 if the order is a delivery)

    • The index of the corresponding delivery order (0 if the order is a pickup)

Program

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

The situation of the pickup-and-delivery problem induces three new constraints:

  • Each pickup/delivery pair must be assigned to the same tour

  • Within this tour, the pickup cannot be visited after the delivery (we cannot deliver an item that has not been picked up yet)

  • The weight of the truck will vary along the tour (increasing upon pickups and decreasing upon deliveries) and the capacity constraint must be satisfied at any time during the tour.

The first requirement for a pickup/delivery pair (a,b) means that a and b must belong to the same list. Literally it means that find(customersSequences, a) == find(customersSequences, b).

As for the ordering of the pair within the list, it is related to the position of a and b within the list, which can be obtained with the indexOf operator: constraint indexOf(pickupList, a) < indexOf(deliveryList, b), with pickupList <- find(customersSequences, a) and deliveryList <- find(customersSequences, b). Note that the first pickup/delivery constraint above makes that pickupList = deliveryList.

As mentioned above, the last specificity is that the capacity constraint must now be checked after each visit. To do this, we need to define the weight after a visit as the weight before the visit plus the weight of the current item (positive if loaded, negative if delivered).

The objectives are the same as for the CVRTW problem: we minimize the total lateness, the number of trucks used, and the total distance traveled by all the trucks.

Execution:
localsolver pdptw.lsp inFileName=instances/lc101.txt [lsTimeLimit=] [solFileName=]
use io;

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

    if (inFileName == nil) throw usage;

    readInputPdptw();

    // Compute distance matrix
    computeDistanceMatrix();

    if (nbMaxTrucks != nil) {
        nbTrucks = nbMaxTrucks;
    }

}

/* 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]);

    //Pickups and deliveries
    for[i in 0...nbCustomers : pickupIndex[i] == -1] {
        pickupListIndex <- find(customersSequences, i);
        deliveryListIndex <- find(customersSequences, deliveryIndex[i]);
        constraint pickupListIndex == deliveryListIndex;
        pickupList <- customersSequences[pickupListIndex];
        deliveryList <- customersSequences[deliveryListIndex];
        constraint indexOf(pickupList, i) < indexOf(deliveryList, deliveryIndex[i]);
    }

    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
        // at any point in the sequence
        routeQuantity[k] <- array(0...c, (i, prev) => prev + demands[sequence[i]], 0);
        constraint and(0...c, i => routeQuantity[k][i] <= truckCapacity);

        endTime[k] <- array(0...c, (i, prev) => max(earliestStart[sequence[i]],
                i == 0 ? distanceDepot[sequence[0]] :
                prev + distanceMatrix[sequence[i - 1]][sequence[i]])
                + serviceTime[sequence[i]], 0);

        homeLateness[k] <- truckUsed[k] ?
               max(0, endTime[k][c - 1] + distanceDepot[sequence[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);

        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]);

    nbTrucksUsed <- sum[k in 0...nbTrucks](truckUsed[k]);

    // Total distance traveled
    totalDistance <- round(100 * sum[k in 0...nbTrucks](routeDistances[k])) / 100;


    minimize totalLateness;
    if (nbMaxTrucks == nil) minimize nbTrucksUsed;
    minimize totalDistance;
}

/* Parametrize 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 readInputPdptw() {
    local inFile = io.openRead(inFileName);

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

    // 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) != 9) 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];
        pickupIndex[i] = line[7].toInt() - 1;
        deliveryIndex[i] = line[8].toInt() - 1;
        i = i + 1;
    }
    nbCustomers = i;

    inFile.close();
}

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

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

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));
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python pdptw.py instances\lc101.txt
Execution (Linux)
export PYTHONPATH=/opt/localsolver_13_0/bin/python
python pdptw.py instances/lc101.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, sol_file):
    #
    # Read instance data
    #
    nb_customers, nb_trucks, truck_capacity, dist_matrix_data, dist_depot_data, \
        demands_data, service_time_data, earliest_start_data, latest_end_data, \
        pick_up_index, delivery_index, max_horizon = read_input_pdptw(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 "at" operators
        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(dist_matrix_data)
        dist_depot = model.array(dist_depot_data)

        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)

        # Pickups and deliveries
        customers_sequences_array = model.array(customers_sequences)
        for i in range(nb_customers):
            if pick_up_index[i] == -1:
                pick_up_list_index = model.find(customers_sequences_array, i)
                delivery_list_index = model.find(customers_sequences_array, delivery_index[i])
                model.constraint(pick_up_list_index == delivery_list_index)
                pick_up_list = model.at(customers_sequences_array, pick_up_list_index)
                delivery_list = model.at(customers_sequences_array, delivery_list_index)
                model.constraint(model.index(pick_up_list, i) < model.index(delivery_list, delivery_index[i]))

        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 at any
            # point in the sequence
            demand_lambda = model.lambda_function(
                lambda i, prev: prev + demands[sequence[i]])
            route_quantity = model.array(model.range(0, c), demand_lambda, 0)

            quantity_lambda = model.lambda_function(
                lambda i: route_quantity[i] <= truck_capacity)
            model.constraint(model.and_(model.range(0, c), quantity_lambda))

            # 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
            end_lambda = model.lambda_function(
                lambda i, prev:
                    model.max(
                        earliest[sequence[i]],
                        model.iif(
                            i == 0,
                            dist_depot[sequence[0]],
                            prev + model.at(dist_matrix, sequence[i - 1], sequence[i])))
                    + service_time[sequence[i]])

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

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

            # Completing visit after latest_end
            late_selector = 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_selector)

        # Total lateness (must be 0 for the solution to be valid)
        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 sol_file is not None:
            with open(sol_file, 'w') as f:
                f.write("%d %.2f\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. +2 is to put it back in
                    # 2...nbCustomers+2 as in the data files (1 being the depot)
                    for customer in customers_sequences[k].value:
                        f.write("%d " % (customer + 1))
                    f.write("\n")


# The input files follow the "Li & Lim" format
def read_input_pdptw(filename):
    file_it = iter(read_elem(filename))

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

    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))

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

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

    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))
        pick = int(next(file_it))
        delivery = 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)
        pick_up_index.append(pick - 1)
        delivery_index.append(delivery - 1)

    nb_customers = i + 1

    distance_matrix = compute_distance_matrix(customers_x, customers_y)
    distance_depots = compute_distance_depots(depot_x, depot_y, customers_x, customers_y)

    return nb_customers, nb_trucks, truck_capacity, distance_matrix, distance_depots, \
            demands, service_time, earliest_start, latest_end, pick_up_index, \
            delivery_index, max_horizon


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


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


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


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

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

    main(instance_file, str_time_limit, sol_file)
Compilation / Execution (Windows)
cl /EHsc pdptw.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver130.lib
pdptw instances\lc101.txt
Compilation / Execution (Linux)
g++ pdptw.cpp -I/opt/localsolver_13_0/include -llocalsolver130 -lpthread -o pdptw
./pdptw instances/lc101.txt
#include "localsolver.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace localsolver;
using namespace std;

class Pdptw {
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;

    // Index for pickup for each node
    vector<int> pickUpIndex;

    // Index for delivery for each node
    vector<int> deliveryIndex;

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

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

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

    Pdptw() {}

    // Read instance data
    void readInstance(const string& fileName) { readInputPdptw(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 "at" operators
        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 distDepot = model.array(distDepotData.begin(), distDepotData.end());

        trucksUsed.resize(nbTrucks);
        vector<LSExpression> distRoutes(nbTrucks), endTime(nbTrucks), homeLateness(nbTrucks), lateness(nbTrucks);

        // Pickups and deliveries
        LSExpression customersSequencesArray = model.array(customersSequences.begin(), customersSequences.end());
        for (int i = 0; i < nbCustomers; ++i) {
            if (pickUpIndex[i] == -1) {
                LSExpression pickUpListIndex = model.find(customersSequencesArray, i);
                LSExpression deliveryListIndex = model.find(customersSequencesArray, deliveryIndex[i]);
                model.constraint(pickUpListIndex == deliveryListIndex);
                LSExpression pickupList = model.at(customersSequencesArray, pickUpListIndex);
                LSExpression deliveryList = model.at(customersSequencesArray, deliveryListIndex);
                model.constraint(model.indexOf(pickupList, i) < model.indexOf(deliveryList, deliveryIndex[i]));
            }
        }

        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 at any point in the sequence
            LSExpression demandLambda = model.createLambdaFunction(
                [&](LSExpression i, LSExpression prev) { return prev + demands[sequence[i]]; });
            LSExpression routeQuantity = model.array(model.range(0, c), demandLambda, 0);

            LSExpression quantityLambda =
                model.createLambdaFunction([&](LSExpression i) { return routeQuantity[i] <= truckCapacity; });
            model.constraint(model.and_(model.range(0, c), quantityLambda));

            // 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
            LSExpression endLambda = model.createLambdaFunction([&](LSExpression i, LSExpression prev) {
                return model.max(earliest[sequence[i]],
                                 model.iif(i == 0, distDepot[sequence[0]],
                                           prev + model.at(distMatrix, sequence[i - 1], sequence[i]))) +
                       serviceTime[sequence[i]];
            });

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

            // Arriving home after max_horizon
            homeLateness[k] =
                model.iif(trucksUsed[k], model.max(0, endTime[k][c - 1] + distDepot[sequence[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
        totalDistance = model.round(100 * model.sum(distRoutes.begin(), distRoutes.end())) / 100;

        // Objective: minimize the number of trucks used, then minimize 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. +2 is to put it back in 2...nbCustomers+2
            // as in the data files (1 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 "Li & Lim" format
    void readInputPdptw(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }

        string str;
        long dump;
        int depotX, depotY;
        vector<int> customersX;
        vector<int> customersY;

        infile >> nbTrucks;
        infile >> truckCapacity;
        infile >> dump;
        infile >> dump;
        infile >> depotX;
        infile >> depotY;
        infile >> dump;
        infile >> dump;
        infile >> maxHorizon;
        infile >> dump;
        infile >> dump;
        infile >> dump;

        while (infile >> dump) {
            int cx, cy, demand, ready, due, service, pick, delivery;
            infile >> cx;
            infile >> cy;
            infile >> demand;
            infile >> ready;
            infile >> due;
            infile >> service;
            infile >> pick;
            infile >> delivery;

            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);
            pickUpIndex.push_back(pick - 1);
            deliveryIndex.push_back(delivery - 1);
        }

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

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

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

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

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

    try {
        Pdptw 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 Pdptw.cs /reference:localsolvernet.dll
Pdptw instances\lc101.txt
using System;
using System.IO;
using System.Collections.Generic;
using localsolver;

public class Pdptw : 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;

    // Index for pick up for each node
    List<int> pickUpIndex;

    // Index for delivery for each node
    List<int> deliveryIndex;

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

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    LSExpression[] customersSequences;

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

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

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

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

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

    // 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 Pdptw()
    {
        localsolver = new LocalSolver();
    }

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputPdptw(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];
        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] = 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 "at" operators
        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 customersSequencesArray = model.Array(customersSequences);

        for (int i = 0; i < nbCustomers; ++i)
        {
            if (pickUpIndex[i] == -1)
            {
                LSExpression pickUpListIndex = model.Find(customersSequencesArray, i);
                LSExpression deliveryListIndex = model.Find(customersSequencesArray, deliveryIndex[i]);
                model.Constraint(pickUpListIndex == deliveryListIndex);

                LSExpression pickUpList = customersSequencesArray[pickUpListIndex];
                LSExpression deliveryList = customersSequencesArray[deliveryListIndex];
                model.Constraint(
                    model.IndexOf(pickUpList, i) < model.IndexOf(deliveryList, deliveryIndex[i])
                );
            }
        }


        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 at any point in the sequence
            LSExpression demandLambda = model.LambdaFunction(
                (i, prev) => prev + demands[sequence[i]]
            );
            LSExpression routeQuantity = model.Array(model.Range(0, c), demandLambda, 0);

            LSExpression quantityLambda = model.LambdaFunction(
                i => routeQuantity[i] <= truckCapacity
            );
            model.Constraint(model.And(model.Range(0, c), quantityLambda));

            // 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
            LSExpression endLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        earliest[sequence[i]],
                        model.If(
                            i == 0,
                            distDepot[sequence[0]],
                            prev + distMatrix[sequence[i - 1], sequence[i]]
                        )
                    ) + serviceTime[sequence[i]]
            );

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

            // Arriving home after max_horizon
            homeLateness[k] = model.If(
                trucksUsed[k],
                model.Max(0, endTime[k][c - 1] + distDepot[sequence[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);
        }

        totalLateness = model.Sum(lateness);
        nbTrucksUsed = model.Sum(trucksUsed);
        totalDistance = model.Round(100 * model.Sum(distRoutes)) / 100;

        // Objective: minimize the number of trucks used, then minimize 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. +2 is to put it back in 2...nbCustomers+2
                // as in the data files (1 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: Pdptw 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 (Pdptw model = new Pdptw())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

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

    // The input files follow the "Li & Lim" format
    private void ReadInputPdptw(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;

            splitted = SplitInput(input);

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

            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>();
            pickUpIndex = new List<int>();
            deliveryIndex = new List<int>();

            while (!input.EndOfStream)
            {
                splitted = SplitInput(input);

                if (splitted.Length < 9)
                    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]);
                pickUpIndex.Add(int.Parse(splitted[7]) - 1);
                deliveryIndex.Add(int.Parse(splitted[8]) - 1);

                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
    private void ComputeDistanceMatrix(
        int depotX,
        int depotY,
        List<int> customersX,
        List<int> customersY
    )
    {
        distMatrixData = new double[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new double[nbCustomers];

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

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

    private double ComputeDist(int xi, int xj, int yi, int yj)
    {
        return Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
    }
}
Compilation / Execution (Windows)
javac Pdptw.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Pdptw instances\lc101.txt
Compilation / Execution (Linux)
javac Pdptw.java -cp /opt/localsolver_13_0/bin/localsolver.jar
java -cp /opt/localsolver_13_0/bin/localsolver.jar:. Pdptw instances/lc101.txt
import java.util.*;
import java.io.*;
import localsolver.*;

public class Pdptw {
    // 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;

    // Index for pick up for each node
    List<Integer> pickUpIndex;

    // Index for delivery for each node
    List<Integer> deliveryIndex;

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

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

    // 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 Pdptw() {
        localsolver = new LocalSolver();
    }

    // Read instance data
    private void readInstance(String fileName) throws IOException {
        readInputPdptw(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 "at" operators
        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);

        // Pickups and deliveries
        LSExpression customersSequencesArray = m.array(customersSequences);
        for (int i = 0; i < nbCustomers; ++i) {
            if (pickUpIndex.get(i) == -1) {
                LSExpression pickUpListIndex = m.find(customersSequencesArray, i);
                LSExpression deliveryListIndex = m.find(customersSequencesArray, deliveryIndex.get(i));
                m.constraint(m.eq(pickUpListIndex, deliveryListIndex));
                LSExpression pickUpList = m.at(customersSequencesArray, pickUpListIndex);
                LSExpression deliveryList = m.at(customersSequencesArray, deliveryListIndex);
                m.constraint(m.leq(m.indexOf(pickUpList, i), m.indexOf(deliveryList, deliveryIndex.get(i))));
            }
        }


        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 at any point in the sequence
            LSExpression demandLambda = m.lambdaFunction((i, prev) -> m.sum(prev, m.at(demands, m.at(sequence, i))));
            LSExpression routeQuantity = m.array(m.range(0, c), demandLambda, 0);

            LSExpression quantityLambda = m.lambdaFunction(i -> m.leq(m.at(routeQuantity, i), truckCapacity));
            m.constraint(m.and(m.range(0, c), quantityLambda));

            // 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
            LSExpression endLambda = 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(distDepot, m.at(sequence, 0)),
                        m.sum(prev, m.at(distMatrix, m.at(sequence, m.sub(i, 1)), m.at(sequence, i)))))),
                m.at(serviceTime, m.at(sequence, i))));

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

            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(distDepot, m.at(sequence, 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. +2 is to put it back in 2...nbCustomers+2
                // as in the data files (1 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 "Li & Lim" format
    private void readInputPdptw(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            nbTrucks = input.nextInt();
            truckCapacity = input.nextInt();
            input.nextInt();

            input.nextInt();
            int depotX = input.nextInt();
            int depotY = input.nextInt();
            input.nextInt();
            input.nextInt();
            maxHorizon = input.nextInt();
            input.nextInt();
            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>();
            pickUpIndex = new ArrayList<Integer>();
            deliveryIndex = 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();
                int pick = input.nextInt();
                int delivery = 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);
                pickUpIndex.add(pick - 1);
                deliveryIndex.add(delivery - 1);
            }

            nbCustomers = customersX.size();

            computeDistanceMatrix(depotX, depotY, customersX, customersY);

        }
    }

    // Compute the distance matrix
    private void computeDistanceMatrix(int depotX, int depotY, List<Integer> customersX, List<Integer> customersY) {
        distMatrixData = new double[nbCustomers][nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 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;
            }
        }

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

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

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

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

            Pdptw model = new Pdptw();
            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);
        }
    }
}