View all code templates

Capacitated Vehicle Routing Problem with Time Windows and Regular Breaks

Routing

Problem

In the Capacitated Vehicle Routing Problem with Time Windows and Regular Breaks, a fleet of delivery vehicles with uniform capacity must service customers. The customers have known opening hours and demand for a single commodity. The vehicles start and end their routes at a common depot, and must take breaks after driving for a certain fixed amount of time. Each customer must be served by exactly one vehicle within its opening hours, and the total demand served by each vehicle must not exceed its capacity. The objectives are to minimize the fleet size and the total traveled distance.

Principles learned

Data

The Vehicle Routing Problem with Time Windows and Regular Breaks instances we provide 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, for 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

Model

The Hexaly model for the Capacitated Vehicle Routing Problem with Time Windows and Regular Breaks extends the CVRPTW model. We refer the reader to this model for the routing and time-windows aspects of the problem. The specific feature of this variant lies in the introduction of mandatory regular breaks, with fixed duration, that must be taken at least every breakFrequency time units.

To model this, we introduce integer decision variables representing the time between successive breaks for each truck. By the break frequency as an upper bound for these decisions, we ensure that breaks are evenly distributed over the planning horizon. The actual break times are then derived cumulatively, taking into account the fixed break duration.

The integration of breaks into route timing adheres to the following principles. Whenever a break occurs during a travel segment, a waiting period, or a service, its fixed duration is added to the current time, thereby shifting all subsequent events along the route.

Finally, the objective remains unchanged from the CVRPTW: we minimize the total lateness, the number of trucks used, and the total traveled distance in lexicographic order.

Download
Execution
hexaly cvrptwrb.hxm inFileName=instances/C101.25.txt [solFileName=] [hxTimeLimit=]
// Copyright (c) Hexaly. Permission is hereby granted to use, copy,
// and modify this code for applications developed with Hexaly.
use io;

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

    if (inFileName == nil) throw usage;

    readInputCvrptwrb();

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

    // A break of 15 minutes every 4 hours
    BREAKFREQUENCY = 60 * 4; //In minutes
    BREAKDURATION = 15;   //In minutes
    nbBreaks = ceil(maxHorizon / BREAKFREQUENCY) + 1;

    // Time between the end of one break and the start of the next
    breaksGaps[k in 0...nbTrucks][p in 0...nbBreaks] <- int(1, BREAKFREQUENCY);
    // Starting time of each break
    breaksStartTimes[k in 0...nbTrucks][b in 0...nbBreaks] <-
        sum[breakIdx in 0...b + 1](breaksGaps[k][breakIdx]) + BREAKDURATION * b;

    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;
        // Breaks must cover the entire horizon
        constraint breaksStartTimes[k][nbBreaks-1] >= maxHorizon + 1;

        // End of each visit
        endTime[k] <- array(0...c, (i, prev) =>
            waitingAndServiceEnd(k, sequence[i], travelEnd(k, i, prev)), 0);

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

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

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

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

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

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

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

/* 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 {trucknumber}: the customers visited [starting and ending service time] | B(starting and ending times) */
function output() {

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

    outfile.println("Instance: ", inFileName);
    outfile.println("Number of trucks: ", nbTrucksUsed.value, " Total distance: ", totalDistance.value, " Max horizon: ", maxHorizon,
                    "\nBreak frequency: ", BREAKFREQUENCY, " Break duration: ", BREAKDURATION, " Working time: ", serviceTime[1]);
    outfile.println("Legend: Client[Start,end] B=Break(Start,end)\n");

    for [k in 0...nbTrucks] {
        if (truckUsed[k].value != 1) continue;
        outfile.print(k, ": ");

        prevEndTime = 0;
        customerOrder = 0;
        customerEndTime = 0;
        customerStartTime = 0;
        for [customer in customersSequences[k].value] {
            customerEndTime = round(endTime[k].value[customerOrder]);
            customerStartTime = customerEndTime - serviceTime[customer];

            // Insert breaks
            for[breakIdx in breaksStartTimes[k]]{
                if (breakIdx.value >= prevEndTime && breakIdx.value <= customerEndTime){
                    endBreak = breakIdx.value + BREAKDURATION;
                    outfile.print("B(", breakIdx.value,", ", endBreak, ") ");
                }
            }
            // # 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:
            outfile.print(customer + 1, "[", customerStartTime ,", " ,customerEndTime, "] ");

            prevEndTime = customerEndTime;
            customerOrder += 1;
        }

        // Insert break if needed before returning to depot
        depotArrivingTime = prevEndTime + distanceDepot[customersSequences[k].value[customerOrder - 1]];
        for[breakIdx in breaksStartTimes[k]]{
            if (breakIdx.value >= prevEndTime && breakIdx.value <= depotArrivingTime){
                endBreak = breakIdx.value + BREAKDURATION;
                outfile.print("B(", breakIdx.value,", ", endBreak, ") ");
                depotArrivingTime += BREAKDURATION;
            }
        }
        outfile.print("| ");
        for[breakIdx in breaksStartTimes[k]]{
        if (breakIdx.value > depotArrivingTime){
            outfile.print("B(", breakIdx.value, ")");
            }
        }
        outfile.print("\n");
    }
}

function readInputCvrptwrb() {
    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));
}

    /* Sub functions for modelling  */
function nextAvailableTime(customer, t) {
    return max(t, earliestStart[customer]);
}

function needsBreak(breakStart, start, end) {
    return start <= breakStart && end > breakStart;
}

    /* Next 3 functions calculate the different times, taking breaks into account. */
function travelEnd(vehicle, i, time) {
    //Compute travel end time considering breaks
    local sequence <- customersSequences[vehicle];
    local travelDuration <- i == 0
            ? distanceDepot[sequence[0]]
            : distanceMatrix[sequence[i - 1]][sequence[i]];
    local travelEnd <- time + travelDuration;
    local endWithBreaks <- travelEnd;
    for [p in 0...nbBreaks] {
        endWithBreaks <- needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks)
                                    ? endWithBreaks + BREAKDURATION
                                    : endWithBreaks;
    }
    return endWithBreaks;
}

function waitingAndServiceEnd(vehicle, customer, time) {
    //Compute waiting and service end time considering breaks
    local nextStartWithoutBreak <- nextAvailableTime(customer, time);
    local endWithoutBreak <- nextStartWithoutBreak + serviceTime[customer];
    local endWithBreaks <- endWithoutBreak;
    for [p in 0...nbBreaks] {
        endWithBreaks <- needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks)
                                    ? nextAvailableTime(customer, breaksStartTimes[vehicle][p] + BREAKDURATION)
                                            + serviceTime[customer]
                                    : endWithBreaks;
    }
    return endWithBreaks;
}
function returningHomeTime(vehicle, customer, time) {
    //Compute returning home time considering breaks
    local endWithoutBreak <- time + distanceDepot[customer];
    local endWithBreaks <- endWithoutBreak;
    for [p in 0...nbBreaks] {
        endWithBreaks <- needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks)
                                    ? endWithBreaks + BREAKDURATION
                                    : endWithBreaks;
    }
    return endWithBreaks;
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python cvrptwrb.py instances\C101.25.txt
Execution (Linux)
export PYTHONPATH=/opt/hexaly_15_0/bin/python
python cvrptwrb.py instances/C101.25.txt
# Copyright (c) Hexaly. Permission is hereby granted to use, copy,
# and modify this code for applications developed with Hexaly.
import hexaly.optimizer
import sys
import math


# Breaks parameters
# A break of 15 minutes every 4 hours
BREAKFREQUENCY = 60*4 # In minutes
BREAKDURATION = 15   # In minutes

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_cvrptwrb(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
        trucks_used = [None] * nb_trucks

        # Number of breaks
        nb_breaks = int(math.ceil(max_horizon / BREAKFREQUENCY) + 1)

        # Time between the end of one break and the start of the next
        breaks_gaps = [[model.int(1, BREAKFREQUENCY) for _ in range(nb_breaks)] for _ in range(nb_trucks)]

        # Starting time of each break
        breaks_start_times = [None] * nb_trucks
        for k in range(0, nb_trucks):
            breaks_start_times_truck = [None] * nb_breaks
            for b in range(0, nb_breaks):
                breaks_start_times_truck[b] = model.sum(breaks_gaps[k][breakIdx] for breakIdx in range(b+1))
                + BREAKDURATION * b
                breaks_start_times[k] = breaks_start_times_truck

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

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

            # 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)
                       # Breaks must cover the entire horizon
            model.constraint(model.geq(breaks_start_times[k][nb_breaks-1], max_horizon + 1))

            # End of each visit
            end_time_lambda = model.lambda_function(
                lambda i, prev:
                    waiting_and_service_end(k, sequence[i],
                                        travel_end(k, i, prev, customers_sequences,model, dist_depot, dist_matrix, nb_breaks, breaks_start_times),
                                        model, earliest, service_time, nb_breaks, breaks_start_times))

            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, returning_home_time(k, sequence[c-1],end_time[k][c-1], model, dist_depot,
                                                    nb_breaks, breaks_start_times) - 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 number of trucks used
        nb_trucks_used = model.sum(trucks_used)

        # 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 {trucknumber}: the customers visited [starting and ending service time] | B(starting and ending times)
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("Instance: " + instance_file + "\n")
                f.write("Number of trucks: " + str(nb_trucks_used.value) + " Total distance: " + str(total_distance.value) + " Max horizon: " + str(max_horizon) + "\n")
                f.write("Break frequency: " + str(BREAKFREQUENCY) + " Break duration: " + str(BREAKDURATION) + " Working time: " + str(service_time_data[1]) + "\n")
                f.write("Legend: Client[Start,end] B=Break(Start,end)\n\n")
                for k in range(nb_trucks):
                    if trucks_used[k].value != 1:
                        continue
                    f.write(str(k) + ": ")

                    prev_end_time = 0
                    customer_order = 0
                    customer_end_time = 0
                    customer_start_time = 0
                    for customer in customers_sequences[k].value:
                        customer_end_time = round(end_time[k].value[customer_order])
                        customer_start_time = customer_end_time - service_time_data[customer]

                        # Insert breaks
                        for breakIdx in breaks_start_times[k]:
                            if (breakIdx.value >= prev_end_time and breakIdx.value <= customer_end_time):
                                end_break = breakIdx.value + BREAKDURATION
                                f.write("B(" + str(breakIdx.value) + "," + str(end_break) + ") ")
                        # 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)
                        f.write(str(customer + 1) + "[" + str(customer_start_time) + "," + str(customer_end_time) + "] ")

                        prev_end_time = customer_end_time
                        customer_order += 1

                    # Insert break if needed before returning to depot
                    depot_arriving_time = prev_end_time + dist_depot_data[customers_sequences[k].value[customer_order - 1]]
                    for breakIdx in breaks_start_times[k]:
                        if (breakIdx.value >= prev_end_time and breakIdx.value <= depot_arriving_time):
                            end_break = breakIdx.value + BREAKDURATION
                            f.write("B(" + str(breakIdx.value) + "," + str(end_break) + ") ")
                            depot_arriving_time += BREAKDURATION

                    f.write("| ")
                    for breakIdx in breaks_start_times[k]:
                        if (breakIdx.value > depot_arriving_time):
                            f.write("B(" + str(breakIdx.value) + ")")
                    f.write("\n")

# The input files follow the "Solomon" format
def read_input_cvrptwrb(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))

    # Sub functions for modelling
def next_available_time(customer, t, model, earliest):
    return model.max(t, earliest[customer])

def needs_break(break_start, start, end, model):
    return model.and_(start <= break_start, end > break_start)

# Next 3 functions compute the different times, taking breaks into account

def travel_end(vehicle, i, time, customers_sequences, model, dist_depot, dist_matrix,
               nb_breaks, breaks_start_times):
    # Compute travel end time
    sequence = customers_sequences[vehicle]
    travel_duration = model.iif(i == 0, dist_depot[sequence[0]], dist_matrix[sequence[i-1]][sequence[i]])
    travel_end = time + travel_duration
    end_with_breaks = travel_end
    for p in range(nb_breaks):
        end_with_breaks = model.iif(needs_break(breaks_start_times[vehicle][p], time, end_with_breaks, model),
                                    end_with_breaks + BREAKDURATION,
                                    end_with_breaks)
    return end_with_breaks

def waiting_and_service_end(vehicle, customer, time, model, earliest, service_time,
                            nb_breaks, breaks_start_times):
    # Compute waiting and service end time
    next_start_without_break = next_available_time(customer, time, model, earliest)
    end_without_break = next_start_without_break + service_time[customer]
    end_with_breaks = end_without_break
    for p in range(nb_breaks):
        end_with_breaks = model.iif(needs_break(breaks_start_times[vehicle][p], time, end_with_breaks, model),
                                    next_available_time(customer, breaks_start_times[vehicle][p] + BREAKDURATION, model, earliest)
                                            + service_time[customer],
                                    end_with_breaks)
    return end_with_breaks

def returning_home_time(vehicle, customer, time, model, dist_depot, nb_breaks, breaks_start_times):
    # Compute returning home time
    end_without_break = time + dist_depot[customer]
    end_with_breaks = end_without_break
    for p in range(nb_breaks):
        end_with_breaks = model.iif(needs_break(breaks_start_times[vehicle][p], time, end_with_breaks, model),
                                    end_with_breaks + BREAKDURATION,
                                    end_with_breaks)
    return end_with_breaks

if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python cvrptwrb.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 cvrptwrb.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly150.lib
cvrptwrb instances\C101.25.txt
Compilation / Execution (Linux)
g++ cvrptwrb.cpp -I/opt/hexaly_15_0/include -lhexaly150 -lpthread -o cvrptwrb
./cvrptwrb instances/C101.25.txt
// Copyright (c) Hexaly. Permission is hereby granted to use, copy,
// and modify this code for applications developed with Hexaly.
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

class Cvrptwrp {
public:

    // Breaks parameters
    // A break of 15 minutes every 4 hours
    int BREAKFREQUENCY = 60*4; //In minutes
    int BREAKDURATION = 15;   //In minutes

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

    // Number of breaks
    int nbBreaks;

    // Decision variables
    vector<HxExpression> customersSequences;

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

    // End time array for each truck
    vector<HxExpression> endTime;

    // Time between the end of one break and the start of the next
    vector<vector<HxExpression>> breaksGaps;

    // Starting time of each break
    vector<vector<HxExpression>> breaksStartTimes;

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

    Cvrptwrp() {}

    /* Read instance data */
    void readInstance(const string& fileName) { readInputCvrptwrp(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);
        endTime.resize(nbTrucks);
        vector<HxExpression> distRoutes(nbTrucks), homeLateness(nbTrucks), lateness(nbTrucks);
        nbBreaks = ceil(maxHorizon / BREAKFREQUENCY) + 1;

        breaksGaps.resize(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k){
            breaksGaps[k].resize(nbBreaks);
            for (int b = 0; b < nbBreaks; ++b){
                breaksGaps[k][b] = model.intVar(1, BREAKFREQUENCY);
            }
        }
        breaksStartTimes.resize(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k){
            breaksStartTimes[k].resize(nbBreaks);
            for (int b = 0; b < nbBreaks; ++b){
                breaksStartTimes[k][b] = model.sum(breaksGaps[k][0]);
                if (b > 0){
                    for (int breakIdx = 1; breakIdx < b + 1; ++breakIdx){
                        breaksStartTimes[k][b].addOperand(breaksGaps[k][breakIdx]);
                    }
                }
                breaksStartTimes[k][b].addOperand(BREAKDURATION * b);
            }
        }

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

            // Breaks must cover the entire horizon
            model.constraint(breaksStartTimes[k][nbBreaks-1] >= maxHorizon + 1);

            // 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 waitingAndServiceEnd(k, sequence[i], travelEnd(k, i, prev, model, distMatrix, distDepot),
                                            model, earliest, serviceTime);
            });

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

            // Arriving home after max horizon
            homeLateness[k] = model.iif(
                    trucksUsed[k],
                    model.max(0, returningHomeTime(k,sequence[c - 1], endTime[k][c - 1], model, distDepot) - 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 {trucknumber}: the customers visited [starting and ending service time] | B(starting and ending times) */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << "Instance: " << fileName << "\n";
        outfile << "Number of trucks: " << nbTrucksUsed.getValue() << " Total distance: " << totalDistance.getDoubleValue() << " Max horizon: " << maxHorizon << endl;
        outfile << "Break frequency: " << BREAKFREQUENCY << " Break duration: " << BREAKDURATION << " Working time: " << serviceTimeData[1] << endl;
        outfile << "Legend: Client[Start,end] B=Break(Start,end)\n" << endl;
        for (int k = 0; k < nbTrucks; ++k) {
            if (trucksUsed[k].getValue() != 1)
                continue;
            outfile << k << ": ";

            int prevEndTime = 0;
            int customerOrder = 0;
            int customerEndTime = 0;
            int customerStartTime = 0;
            HxCollection customersCollection = customersSequences[k].getCollectionValue();
            for (int i = 0; i < customersCollection.count(); ++i) {
                int customer = customersCollection[i];
                customerEndTime = round(endTime[k].getArrayValue().getDoubleValue(customerOrder));
                customerStartTime =  customerEndTime - serviceTimeData[customer];

                // Insert breaks
                for (const HxExpression& breakIdx : breaksStartTimes[k]) {
                    if (breakIdx.getValue() >= prevEndTime && breakIdx.getValue() <= customerEndTime){
                        int endBreak = breakIdx.getValue() + BREAKDURATION;
                        outfile << "B(" << breakIdx.getValue() << ", " << endBreak << ") ";
                    }
                }
                // 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)
                outfile << customer + 1 << "[" << customerStartTime << ", " << customerEndTime << "] ";

                prevEndTime = customerEndTime;
                customerOrder += 1;
            }

            // Insert break if needed before returning to depot
            int depotArrivingTime = prevEndTime + distDepotData[customersCollection[customerOrder - 1]];
            for (const HxExpression& breakIdx : breaksStartTimes[k]) {
                if (breakIdx.getValue() >= prevEndTime && breakIdx.getValue() <= depotArrivingTime){
                    int endBreak = breakIdx.getValue() + BREAKDURATION;
                    outfile << "B(" << breakIdx.getValue() << ", " << endBreak << ") ";
                    depotArrivingTime += BREAKDURATION;
                }
            }
            outfile << "| ";
            for (const HxExpression& breakIdx : breaksStartTimes[k]) {
                if (breakIdx.getValue() > depotArrivingTime){
                    outfile << "B(" << breakIdx.getValue() << ")";
                    }
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the "Solomon" format
    void readInputCvrptwrp(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));
    }

    // Sub functions for modelling
    HxExpression nextAvailableTime(HxExpression  customer, HxExpression  t, HxModel model, HxExpression earliest){
        return model.max(t, earliest[customer]);}

    HxExpression needsBreak(HxExpression  breakStart, HxExpression  start, HxExpression  end, HxModel model){
        return model.and_(start <= breakStart, end > breakStart);
    }

    // Next 3 functions compute the different times, taking breaks into account
    HxExpression travelEnd(int vehicle, HxExpression  i, HxExpression time, HxModel model, HxExpression distMatrix, HxExpression distDepot){
        //Compute travel end time
        HxExpression sequence = customersSequences[vehicle];
        HxExpression travelDuration = model.iif(i == 0, distDepot[sequence[0]], distMatrix[sequence[i-1]][sequence[i]]);
        HxExpression travelEnd = time + travelDuration;
        HxExpression endWithBreaks = travelEnd;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    endWithBreaks + BREAKDURATION,
                                    endWithBreaks);
        }
        return endWithBreaks;
    }

    HxExpression waitingAndServiceEnd(int vehicle, HxExpression customer, HxExpression time, HxModel model, HxExpression earliest, HxExpression serviceTime){
        //Compute waiting and service end time
        HxExpression nextStartWithoutBreak = nextAvailableTime(customer, time, model, earliest);
        HxExpression endWithoutBreak = nextStartWithoutBreak + serviceTime[customer];
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    nextAvailableTime(customer, breaksStartTimes[vehicle][p] + BREAKDURATION, model, earliest)
                                            + serviceTime[customer],
                                    endWithBreaks);
        }
        return endWithBreaks;
    }

    HxExpression returningHomeTime(int vehicle, HxExpression  customer, HxExpression  time, HxModel model, HxExpression distDepot){
        //Compute returning home time
        HxExpression endWithoutBreak = time + distDepot[customer];
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    endWithBreaks + BREAKDURATION,
                                    endWithBreaks);
        }
        return endWithBreaks;
    }

};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: cvrptwrp 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 {
        Cvrptwrp 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 Cvrptwrb.cs /reference:Hexaly.NET.dll
Cvrptwrb instances\C101.25.txt
// Copyright (c) Hexaly. Permission is hereby granted to use, copy,
// and modify this code for applications developed with Hexaly.
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

public class Cvrptwrb : IDisposable
{
    // Breaks parameters
    // A break of 15 minutes every 4 hours
    int BREAKFREQUENCY = 60*4; //In minutes
    int BREAKDURATION = 15;   //In minutes
    // 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;

    // Number of breaks
    int nbBreaks;

    // Decision variables
    HxExpression[] customersSequences;

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

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

    // Time between the end of one break and the start of the next
    HxExpression[][] breaksGaps;

    // Starting time of each break
    HxExpression[][] breaksStartTimes;

    // 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 Cvrptwrb()
    {
        optimizer = new HexalyOptimizer();
    }

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputCvrptwrb(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];
        endTime = new HxExpression[nbTrucks];
        HxExpression[] distRoutes = 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);

        //Add breaks
        nbBreaks = (int) Math.Ceiling( (double) maxHorizon/ (double) BREAKFREQUENCY) + 1;

        breaksGaps = new HxExpression[nbTrucks][];
        for (int k = 0; k <nbTrucks; ++k){
            breaksGaps[k] = new HxExpression[nbBreaks];
            for(int b = 0; b < nbBreaks; ++b){
                breaksGaps[k][b] = model.Int(1,BREAKFREQUENCY);
            }
        }
        breaksStartTimes = new HxExpression[nbTrucks][];
        for (int k = 0; k <nbTrucks; ++k){
            breaksStartTimes[k] = new HxExpression[nbBreaks];
            for(int b = 0; b < nbBreaks; ++b){
                breaksStartTimes[k][b] = model.Sum(breaksGaps[k][0]);
                if (b>0){
                    for (int breakIdx = 1; breakIdx <= b; ++breakIdx){
                        breaksStartTimes[k][b].AddOperand(breaksGaps[k][breakIdx]);
                    }
                }
                breaksStartTimes[k][b].AddOperand(BREAKDURATION * b);
            }
        }

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

            // Breaks must cover the entire horizon
            model.Constraint(breaksStartTimes[k][nbBreaks-1] >= maxHorizon + 1);

            // 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) =>
                    waitingAndServiceEnd(k, sequence[i], travelEnd(k, i, prev, model, distMatrix, distDepot), model, earliest, serviceTime)
            );

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

            // Arriving home after max_horizon
            homeLateness[k] = model.If(
                trucksUsed[k],
                model.Max(0, returningHomeTime(k,sequence[c - 1], endTime[k][c - 1], model, distDepot)  - 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 {trucknumber}: the customers visited [starting and ending service time] | B(starting and ending times) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("Instance: " + fileName);
            output.WriteLine("Number of trucks: " + nbTrucksUsed.GetIntValue() + " Total distance: " + totalDistance.GetDoubleValue() + " Max horizon: " + maxHorizon +
                    "\nBreak frequency: " + BREAKFREQUENCY + " Break duration: " + BREAKDURATION + " Working time: " + serviceTimeData[1]);
            output.WriteLine("Legend: Client[Start,end] B=Break(Start,end)\n");
            for (int k = 0; k < nbTrucks; ++k)
            {
                if (trucksUsed[k].GetValue() != 1)
                    continue;
                output.Write(k + ": ");

                int prevEndTime = 0;
                int customerOrder = 0;
                int customerEndTime = 0;
                int customerStartTime = 0;
                // 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) {
                    int  customer = (int) customersCollection[i];
                    customerEndTime = (int) Math.Round( endTime[k].GetArrayValue().GetDoubleValue(customerOrder));
                    customerStartTime = customerEndTime - serviceTimeData[customer];

                    // Insert breaks
                    foreach (HxExpression breakIdx in breaksStartTimes[k]) {
                        if (breakIdx.GetIntValue() >= prevEndTime && breakIdx.GetIntValue() <= customerEndTime){
                            int endBreak = (int) breakIdx.GetValue() + BREAKDURATION;
                            output.Write("B(" + breakIdx.GetValue() + ", " + endBreak + ") ");
                        }
                    }
                    output.Write((customer + 1) + "[" + customerStartTime + ", " + customerEndTime + "] ");
                    prevEndTime = customerEndTime;
                    customerOrder += 1;
                }

                // Insert break if needed before returning to depot
                int depotArrivingTime = prevEndTime + (int) distDepotData[customersCollection[customerOrder - 1]];
                foreach (HxExpression breakIdx in breaksStartTimes[k]) {
                    if (breakIdx.GetIntValue() >= prevEndTime && breakIdx.GetIntValue() <= depotArrivingTime){
                        int endBreak = (int) breakIdx.GetValue() + BREAKDURATION;
                        output.Write("B(" + breakIdx.GetValue() + ", " + endBreak + ") ");
                        depotArrivingTime += BREAKDURATION;
                    }
                }
                output.Write("| ");
                foreach (HxExpression breakIdx in breaksStartTimes[k]){
                    if (breakIdx.GetValue() > depotArrivingTime) {
                        output.Write("B(" + breakIdx.GetValue() + ")");
                    }
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Cvrptwrb 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 (Cvrptwrb model = new Cvrptwrb())
        {
            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 ReadInputCvrptwrb(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));
    }

    /* Breaks functions */
    HxExpression nextAvailableTime(HxExpression  customer, HxExpression  t, HxModel model, HxExpression earliest){
        return model.Max(t, earliest[customer]);}

    HxExpression needsBreak(HxExpression  breakStart, HxExpression  start, HxExpression  end, HxModel model){
        return model.And(start <= breakStart,end > breakStart);
    }

    // Next 3 functions compute the different times, taking breaks into account
    HxExpression travelEnd(int vehicle, HxExpression  i, HxExpression time, HxModel model, HxExpression distMatrix, HxExpression distDepot){
        //Compute travel end time
        HxExpression sequence = customersSequences[vehicle];
        HxExpression travelDuration = model.If(i == 0, distDepot[sequence[0]],distMatrix[sequence[i-1]][sequence[i]]);
        HxExpression travelEnd = time + travelDuration;
        HxExpression endWithBreaks = travelEnd;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.If(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    endWithBreaks + BREAKDURATION,
                                    endWithBreaks);
       }
        return endWithBreaks;
    }

    HxExpression waitingAndServiceEnd(int vehicle, HxExpression customer, HxExpression time, HxModel model, HxExpression earliest, HxExpression serviceTime){
        //Compute waiting and service end time
        HxExpression nextStartWithoutBreak = nextAvailableTime(customer, time, model, earliest);
        HxExpression endWithoutBreak = nextStartWithoutBreak + serviceTime[customer];
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.If(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    nextAvailableTime(customer, breaksStartTimes[vehicle][p] + BREAKDURATION, model, earliest)
                                            + serviceTime[customer],
                                    endWithBreaks);
       }
       return endWithBreaks;
    }

    HxExpression returningHomeTime(int vehicle, HxExpression  customer, HxExpression  time, HxModel model, HxExpression distDepot){
        //Compute returning home time
        HxExpression endWithoutBreak = time + distDepot[customer];
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = model.If(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, model),
                                    endWithBreaks + BREAKDURATION,
                                    endWithBreaks);
       }
       return endWithBreaks;
    }
}
Compilation / Execution (Windows)
javac Cvrptwrb.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. Cvrptwrb instances\C101.25.txt
Compilation / Execution (Linux)
javac Cvrptwrb.java -cp /opt/hexaly_15_0/bin/hexaly.jar
java -cp /opt/hexaly_15_0/bin/hexaly.jar:. Cvrptwrb instances/C101.25.txt
// Copyright (c) Hexaly. Permission is hereby granted to use, copy,
// and modify this code for applications developed with Hexaly.
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class Cvrptwrb {

    // Break parameters
    // A break of 15 minutes every 4 hours
    private int BREAKFREQUENCY = 60*4; //In minutes
    private int BREAKDURATION = 15;   //In minutes

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

    // Number of breaks
    private int nbBreaks;

    // 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;
    
    // Time between the end of one break and the start of the next
    private HxExpression[][] breaksGaps;

    // Starting time of each break
    private HxExpression[][] breaksStartTimes;

    // 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 Cvrptwrb(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

    /* Read instance data */
    private void readInstance(String fileName) throws IOException {
        readInputCvrptwrb(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);

        nbBreaks = (int) Math.ceil(maxHorizon / BREAKFREQUENCY) + 1;
        breaksGaps = new HxExpression[nbTrucks][nbBreaks];
        for (int k = 0; k < nbTrucks; ++k){
            for(int b = 0; b < nbBreaks; ++b){
                breaksGaps[k][b] = m.intVar(1, BREAKFREQUENCY);
            }
        }
        breaksStartTimes = new HxExpression[nbTrucks][nbBreaks];
        for (int k = 0; k <nbTrucks; ++k){
            for(int b = 0; b < nbBreaks; ++b){
                breaksStartTimes[k][b] = m.sum(breaksGaps[k][0]);
                if (b > 0){
                    for (int breakIdx = 1; breakIdx <= b; ++breakIdx){
                        breaksStartTimes[k][b].addOperand(breaksGaps[k][breakIdx]);
                    }
                }
                breaksStartTimes[k][b].addOperand(BREAKDURATION * b);
            }
        }

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

            // Breaks must cover the entire horizon
            m.constraint(m.geq(breaksStartTimes[k][nbBreaks-1], maxHorizon + 1));

            // 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
            int truck = k;
            HxExpression endTimeLambda = m.lambdaFunction((i, prev) ->
                waitingAndServiceEnd(truck, m.at(sequence, i), travelEnd(truck, i, prev, m, distMatrix, distDepot),
                                        m, earliest, serviceTime));

            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.sub(returningHomeTime(k, m.at(sequence, m.sub(c, 1)), m.at(endTime[k],m.sub(c, 1)), m, distDepot), maxHorizon)),
                0
            );

            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 {trucknumber}: the customers visited [starting and ending service time] | B(starting and ending times)
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println("Instance: " + fileName);
            output.println("Number of trucks: " + nbTrucksUsed.getIntValue() + " Total distance: " + totalDistance.getDoubleValue() + " Max horizon: " + maxHorizon +
                    "\nBreak frequency: " + BREAKFREQUENCY + " Break duration: " + BREAKDURATION + " Working time: " + serviceTimeData.get(1));
            output.println("Legend: Client[Start,end] B=Break(Start,end)\n");
                for (int k = 0; k < nbTrucks; ++k) {
                if (trucksUsed[k].getValue() != 1)
                    continue;
                output.print(k + ": ");

                int prevEndTime = 0;
                int customerOrder = 0;
                int customerEndTime = 0;
                int customerStartTime = 0;
                HxCollection customersCollection = customersSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); ++i) {
                    int  customer = (int) customersCollection.get(i);
                    customerEndTime = Math.round((int) endTime[k].getArrayValue().getDoubleValue(customerOrder));
                    customerStartTime = customerEndTime - serviceTimeData.get(customer);

                    // Insert breaks
                    for (HxExpression breakIdx : breaksStartTimes[k]) {
                        if (breakIdx.getValue() >= prevEndTime && breakIdx.getValue() <= customerEndTime){
                            int endBreak = (int) breakIdx.getValue() + BREAKDURATION;
                            output.print("B(" + breakIdx.getValue() + ", " + endBreak + ") ");
                        }
                    }
                    // 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)
                    output.print((customer + 1) + "[" + customerStartTime + ", " + customerEndTime + "] ");

                    prevEndTime = customerEndTime;
                    customerOrder += 1;
                }

                // Insert break if needed before returning to depot
                int depotArrivingTime = prevEndTime + (int) distDepotData[(int) customersCollection.get(customerOrder - 1)];
                for (HxExpression breakIdx : breaksStartTimes[k]) {
                    if (breakIdx.getValue() >= prevEndTime && breakIdx.getValue() <= depotArrivingTime){
                        int endBreak = (int) breakIdx.getValue() + BREAKDURATION;
                        output.print("B(" + breakIdx.getValue() + ", " + endBreak + ") ");
                        depotArrivingTime += BREAKDURATION;
                    }
                }
                output.print("| ");
                for (HxExpression breakIdx : breaksStartTimes[k]){
                    if (breakIdx.getValue() > depotArrivingTime) {
                        output.print("B(" + breakIdx.getIntValue() + ")");
                    }
                }
                output.print("\n");
            }
        }
    }


    // The input files follow the "Solomon" format
    private void readInputCvrptwrb(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.useLocale(Locale.ROOT);
            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));
    }

    // Sub functions for modelling
    private HxExpression nextAvailableTime(HxExpression  customer, HxExpression  t, HxModel m, HxExpression earliest){
        return m.max(t, m.at(earliest, customer));}

    private HxExpression needsBreak(HxExpression  breakStart, HxExpression  start, HxExpression  end, HxModel m){
        return m.and(m.leq(start, breakStart),m.gt(end, breakStart));
    }

    // Next 3 functions compute the different times, taking breaks into account
    private HxExpression travelEnd(int vehicle, HxExpression  i, HxExpression time, HxModel m,
                                    HxExpression distMatrix, HxExpression distDepot){
        //Compute travel end time
        HxExpression sequence = customersSequences[vehicle];
        HxExpression travelDuration = m.iif(m.eq(i, 0), m.at(distDepot, m.at(sequence, 0)),
                                            m.at(distMatrix, m.at(sequence, m.sub(i,1)), m.at(sequence, i)));
        HxExpression travelEnd = m.sum(time, travelDuration);
        HxExpression endWithBreaks = travelEnd;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = m.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, m),
                                    m.sum(endWithBreaks, BREAKDURATION),
                                    endWithBreaks);
       }
        return endWithBreaks;
    }

    private HxExpression waitingAndServiceEnd(int vehicle, HxExpression customer, HxExpression time, HxModel m,
                                                HxExpression earliest, HxExpression serviceTime){
        //Compute waiting and service end time
        HxExpression nextStartWithoutBreak = nextAvailableTime(customer, time, m, earliest);
        HxExpression endWithoutBreak = m.sum(nextStartWithoutBreak, m.at(serviceTime,customer));
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = m.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, m),
                                    m.sum(nextAvailableTime(customer, m.sum(breaksStartTimes[vehicle][p],BREAKDURATION), m, earliest),
                                                m.at(serviceTime, customer)),
                                    endWithBreaks);
       }
       return endWithBreaks;
    }

    private HxExpression returningHomeTime(int vehicle, HxExpression  customer, HxExpression  time, HxModel m, HxExpression distDepot){
        //Compute returning home time
        HxExpression endWithoutBreak = m.sum(time, m.at(distDepot,customer));
        HxExpression endWithBreaks = endWithoutBreak;
        for(int p = 0; p < nbBreaks; ++p){
            endWithBreaks = m.iif(needsBreak(breaksStartTimes[vehicle][p], time, endWithBreaks, m),
                                    m.sum(endWithBreaks, BREAKDURATION),
                                    endWithBreaks);
       }
       return endWithBreaks;
    }

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

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

Explore templates

Discover the ease of use and performance of Hexaly