Vehicle Routing with Time Windows (CVRPTW)¶

Principles learned¶

  • Add multiple list decision variables

  • Use lambda expressions to define a recursive array

  • Define a sequence of expressions

Problem¶

../_images/vrptw.svg

In the capacitated vehicle routing problem with time-windows (CVRPTW), 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. 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.

Download the example


Data¶

The instances provided come from the Solomon 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

Program¶

The Hexaly model is an extension of the CVRP model. We refer the reader to this model for the routing aspect of the problem. The time-windows are handled as a prioritary objective. In case of early arrival, the truck has to wait for the opening hour of the customer. In case of late arrival, the lateness is measured and the total lateness has to be minimized. The solution is feasible when this cumulated lateness is zero.

In the model the ending times of each truck are defined as a recursive array. The arrival time is the maximum of:

  • the end time of the previous visit plus the traveling time (equal to the distance in these instances). For the first visit it amounts to the traveling time from the depot (starting at t=0)

  • the earliest allowed arrival at this customer

The ending time is merely this arrival time plus the service time for this customer. The arrival to the depot at the end of the tour occurs at the end time of the last visit plus the traveling time from this point to the depot.

Such a recursive definition is introduced with a function (i, prev) => .... It allows defining the ith element of an array in function of i and of the (i-1)th element of the array. See our documentation on this topic for details.

The lateness at each visit is computed from the difference between the ending time of this visit and the latest allowed end for this customer. For the arrival at the depot, we compare the arrival time to the maximum horizon defined for the problem.

Finally we minimize in lexicographic order: the total lateness, the number of trucks used, and the total distance traveled by all the trucks.

Execution:
hexaly cvrptw.hxm inFileName=instances/R101.25.txt [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data. The input files follow the "Solomon" format*/
function input() {
    usage = "Usage: hexaly cvrptw.hxm "
            + "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;

    readInputCvrptw();

    computeDistanceMatrix();
}

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

        // Arriving home after max horizon
        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);

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


/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 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 matrix
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=%HX_HOME%\bin\python
python cvrptw.py instances\R101.25.txt
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python cvrptw.py instances/R101.25.txt
import hexaly.optimizer
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, dist_depot_data, \
        demands_data, service_time_data, earliest_start_data, latest_end_data, \
        max_horizon = read_input_cvrptw(instance_file)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        model = optimizer.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 Hexaly 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(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)

        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
            end_time_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_time_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_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 optimizer
        optimizer.param.time_limit = int(str_time_limit)

        optimizer.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

    # Compute distance matrix
    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, max_horizon


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


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


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


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python cvrptw.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 cvrptw.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
cvrptw instances\R101.25.txt
Compilation / Execution (Linux)
g++ cvrptw.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o cvrptw
./cvrptw instances/R101.25.txt
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

class Cvrptw {
public:
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    vector<HxExpression> customersSequences;

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

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

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    Cvrptw() {}

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

    void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.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 Hexaly arrays to be able to access them with an "at" operator
        HxExpression demands = model.array(demandsData.begin(), demandsData.end());
        HxExpression earliest = model.array(earliestStartData.begin(), earliestStartData.end());
        HxExpression latest = model.array(latestEndData.begin(), latestEndData.end());
        HxExpression serviceTime = model.array(serviceTimeData.begin(), serviceTimeData.end());
        HxExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            distMatrix.addOperand(model.array(distMatrixData[n].begin(), distMatrixData[n].end()));
        }
        HxExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());

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

        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = customersSequences[k];
            HxExpression 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
            HxExpression demandLambda =
                model.createLambdaFunction([&](HxExpression j) { return demands[j]; });
            HxExpression routeQuantity = model.sum(sequence, demandLambda);
            model.constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            HxExpression distLambda = model.createLambdaFunction(
                [&](HxExpression 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
            HxExpression endTimeLambda = model.createLambdaFunction([&](HxExpression i, HxExpression 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), endTimeLambda, 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
            HxExpression lateLambda = model.createLambdaFunction(
                [&](HxExpression 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();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.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)
            HxCollection 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);
        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: cvrptw 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 {
        Cvrptw 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 %HX_HOME%\bin\Hexaly.NET.dll .
csc Cvrptw.cs /reference:Hexaly.NET.dll
Cvrptw instances\R101.25.txt
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

public class Cvrptw : IDisposable
{
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    HxExpression[] customersSequences;

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

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

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    public Cvrptw()
    {
        optimizer = new HexalyOptimizer();
    }

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

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

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

        trucksUsed = new HxExpression[nbTrucks];
        customersSequences = new HxExpression[nbTrucks];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] endTime = new HxExpression[nbTrucks];
        HxExpression[] homeLateness = new HxExpression[nbTrucks];
        HxExpression[] lateness = new HxExpression[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 HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.Array(demandsData);
        HxExpression earliest = model.Array(earliestStartData);
        HxExpression latest = model.Array(latestEndData);
        HxExpression serviceTime = model.Array(serviceTimeData);
        HxExpression distDepot = model.Array(distDepotData);
        HxExpression distMatrix = model.Array(distMatrixData);

        for (int k = 0; k < nbTrucks; ++k)
        {
            HxExpression sequence = customersSequences[k];
            HxExpression 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
            HxExpression demandLambda = model.LambdaFunction(j => demands[j]);
            HxExpression routeQuantity = model.Sum(sequence, demandLambda);
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            HxExpression 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
            HxExpression endTimeLambda = 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), endTimeLambda, 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
            HxExpression 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();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.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)
                HxCollection 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: Cvrptw 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 (Cvrptw model = new Cvrptw())
        {
            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[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
    }

    // The input files follow the "Solomon" format
    private 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
    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 Cvrptw.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. Cvrptw instances\R101.25.txt
Compilation / Execution (Linux)
javac Cvrptw.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. Cvrptw instances/R101.25.txt
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class Cvrptw {
    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

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

    // Number of trucks
    private int nbTrucks;

    // Decision variables
    private HxExpression[] customersSequences;

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

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

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

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

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

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

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

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

    private Cvrptw(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

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

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

        trucksUsed = new HxExpression[nbTrucks];
        customersSequences = new HxExpression[nbTrucks];
        distRoutes = new HxExpression[nbTrucks];
        endTime = new HxExpression[nbTrucks];
        homeLateness = new HxExpression[nbTrucks];
        lateness = new HxExpression[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 HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = m.array(demandsData);
        HxExpression earliest = m.array(earliestStartData);
        HxExpression latest = m.array(latestEndData);
        HxExpression serviceTime = m.array(serviceTimeData);
        HxExpression distDepot = m.array(distDepotData);
        HxExpression distMatrix = m.array(distMatrixData);

        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = customersSequences[k];
            HxExpression 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
            HxExpression demandLambda = m.lambdaFunction(j -> m.at(demands, j));
            HxExpression routeQuantity = m.sum(sequence, demandLambda);
            m.constraint(m.leq(routeQuantity, truckCapacity));

            // Distance traveled by truck k
            HxExpression 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
            HxExpression 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(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), endTimeLambda, 0);

            HxExpression 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
            HxExpression 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();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.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)
                HxCollection 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];
        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 Cvrptw inputFile [outputFile] [timeLimit] [nbTrucks]");
            System.exit(1);
        }

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

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