Time-Dependent Routing Problem with Time Windows (TDCVRPTW)

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, travel_time_data, \
        time_to_matrix_idx_data, dist_depot_data, travel_time_warehouse_data,\
        demands_data, service_time_data, earliest_start_data, latest_end_data, \
        max_horizon = read_input_cvrptw(instance_file)

    with 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)
        travel_time = model.array(travel_time_data)
        time_to_matrix_idx = model.array(time_to_matrix_idx_data)
        dist_depot = model.array(dist_depot_data)
        travel_time_warehouse = model.array(travel_time_warehouse_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 according to the traffic
            end_time_lambda = model.lambda_function(
                lambda i, prev:
                    model.max(
                        earliest[sequence[i]],
                        model.iif(
                            i == 0, model.at(travel_time_warehouse, sequence[0], time_to_matrix_idx[0]),
                            prev + model.at(travel_time, sequence[i - 1], sequence[i],
                                            time_to_matrix_idx[model.round(prev)])))
                    + service_time[sequence[i]])

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

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

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

        # Total lateness
        total_lateness = model.sum(lateness)

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

        # Objective: minimize the number of trucks used, then minimize the distance traveled
        model.minimize(total_lateness)
        model.minimize(nb_trucks_used)
        model.minimize(total_distance)

        model.close()

        # Parameterize the 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

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

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

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


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

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

        for k in range(nb_time_intervals):
            travel_time[i][i][k] = 0

        for j in range(nb_customers):
            dist = compute_dist(customers_x[i], customers_x[j],
                                customers_y[i], customers_y[j])
            distance_matrix[i][j] = dist
            distance_matrix[j][i] = dist

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

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


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

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


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


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


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

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

    main(instance_file, str_time_limit, output_file)

#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

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

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

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

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

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

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

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

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

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    vector<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;

    Tdcvrptw() {}

    /* 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 travelTime = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            HxExpression timeMatrix = model.array();
            for (int m = 0; m < nbCustomers; ++m) {
                timeMatrix.addOperand(model.array(travelTimeData[n][m].begin(), travelTimeData[n][m].end()));
            }
            travelTime.addOperand(timeMatrix);
        }
        HxExpression timeToMatrixIdx = model.array(timeToMatrixIdxData.begin(), timeToMatrixIdxData.end());
        HxExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());
        HxExpression travelTimeWarehouse = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            travelTimeWarehouse.addOperand(
                model.array(travelTimeWarehouseData[n].begin(), travelTimeWarehouseData[n].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 according to the traffic
            HxExpression endTimeLambda = model.createLambdaFunction([&](HxExpression i, HxExpression prev) {
                return model.max(earliest[sequence[i]],
                                 model.iif(i == 0, model.at(travelTimeWarehouse, sequence[0], timeToMatrixIdx[0]),
                                           prev + model.at(travelTime, sequence[i - 1], sequence[i],
                                                           timeToMatrixIdx[model.round(prev)]))) +
                       serviceTime[sequence[i]];
            });

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

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

            // Completing visit after latest end
            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();

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

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

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

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

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

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

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

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

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

    try {
        Tdcvrptw model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}

using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

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

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

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

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

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

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

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

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

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    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 Tdcvrptw()
    {
        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);
        HxExpression travelTime = model.Array(travelTimeData);
        HxExpression timeToMatrixIdx = model.Array(timeToMatrixIdxData);
        HxExpression travelTimeWarehouse = model.Array(travelTimeWarehouseData);

        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 according to the traffic
            HxExpression endTimeLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        earliest[sequence[i]],
                        model.If(
                            i == 0,
                            travelTimeWarehouse[sequence[0], timeToMatrixIdx[0]],
                            prev
                                + model.At(
                                    travelTime,
                                    sequence[i - 1],
                                    sequence[i],
                                    timeToMatrixIdx[model.Round(prev)]
                                )
                        )
                    ) + serviceTime[sequence[i]]
            );

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

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

            // Completing visit after latest_end
            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();

        // Parameterize 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: Tdcvrptw inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

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

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

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

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

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

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

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

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

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

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

            nbCustomers = customersX.Count;

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

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

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

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

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

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

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

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

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

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

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

import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

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

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

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

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

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

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

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

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

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

    // Number of trucks
    private int nbTrucks;

    // Decision variables
    private 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 Tdcvrptw(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);
        HxExpression travelTime = m.array(travelTimeData);
        HxExpression timeToMatrixIdx = m.array(timeToMatrixIdxData);
        HxExpression travelTimeWarehouse = m.array(travelTimeWarehouseData);

        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 according to the traffic
            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(travelTimeWarehouse, m.at(sequence, 0), m.at(timeToMatrixIdx, 0)),
                                    m.sum(prev, m.at(travelTime, m.at(sequence, m.sub(i, 1)), m.at(sequence, i),
                                            m.at(timeToMatrixIdx, m.round(prev))))))),
                    m.at(serviceTime, m.at(sequence, i))));

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

            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(m.at(travelTimeWarehouse, m.at(sequence, m.sub(c, 1))),
                                            m.at(timeToMatrixIdx, m.round(m.at(theEnd, m.sub(c, 1))))), maxHorizon))),
                    0);

            // Completing visit after latest_end
            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();

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

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

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

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

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

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

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

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

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

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