Dial-A-Ride Problem (DARP)

import hexaly.optimizer
import sys
import json

def read_data(filename):
    with open(filename) as f:
        return json.load(f)
    
def read_input_darp(instance_file):
    instance = read_data(instance_file)

    nb_clients = instance["nbClients"]
    nb_nodes = instance["nbNodes"]
    nb_vehicles = instance["nbVehicles"]
    depot_tw_end = instance["depot"]["twEnd"]
    capacity = instance["capacity"]
    scale = instance["scale"]

    quantities = [-1 for i in range(2 * nb_clients)]
    distances = instance["distanceMatrix"]
    starts = [-1.0 for i in range(2 * nb_clients)]
    ends = [-1.0 for i in range(2 * nb_clients)]
    loading_times = [-1.0 for i in range(2 * nb_clients)]
    max_travel_times = [-1.0 for i in range(2 * nb_clients)]
    for k in range(nb_clients):
        quantities[k] = instance["clients"][k]["nbClients"]
        quantities[k+nb_clients] = -instance["clients"][k]["nbClients"]

        starts[k] = instance["clients"][k]["pickup"]["start"]
        ends[k] = instance["clients"][k]["pickup"]["end"]
        
        starts[k+nb_clients] = instance["clients"][k]["delivery"]["start"]
        ends[k+nb_clients] = instance["clients"][k]["delivery"]["end"]

        loading_times[k] = instance["clients"][k]["pickup"]["loadingTime"]
        loading_times[k+nb_clients] = instance["clients"][k]["delivery"]["loadingTime"]

        max_travel_times[k] = instance["clients"][k]["pickup"]["maxTravelTime"]
        max_travel_times[k+nb_clients] = instance["clients"][k]["delivery"]["maxTravelTime"]
        
    factor = 1.0 / (scale * instance["speed"])

    distance_warehouse = [-1.0 for i in range(nb_nodes)]
    time_warehouse = [-1.0 for i in range(nb_nodes)]
    distance_matrix = [[-1.0 for i in range(nb_nodes)] for j in range(nb_nodes)]
    time_matrix = [[-1.0 for i in range(nb_nodes)] for j in range(nb_nodes)]
    for i in range(nb_nodes):
        distance_warehouse[i] = distances[0][i+1]
        time_warehouse[i] = distance_warehouse[i] * factor
        for j in range(nb_nodes):
            distance_matrix[i][j] = distances[i+1][j+1]
            time_matrix[i][j] = distance_matrix[i][j] * factor

    return nb_clients, nb_nodes, nb_vehicles, depot_tw_end, capacity, scale, quantities, \
        starts, ends, loading_times, max_travel_times, distance_warehouse, time_warehouse, \
        distance_matrix, time_matrix

def main(instance_file, str_time_limit, sol_file):

    nb_clients, nb_nodes, nb_vehicles, depot_tw_end, capacity, scale, quantities_data, \
        starts_data, ends_data, loading_times_data, max_travel_times, distance_warehouse_data, \
        time_warehouse_data, distance_matrix_data, time_matrix_data = read_input_darp(instance_file)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        model = optimizer.model

        # routes[k] represents the nodes visited by vehicle k
        routes = [model.list(nb_nodes) for k in range(nb_vehicles)]
        depot_starts = [model.float(0, depot_tw_end) for k in range(nb_vehicles)]
        # Each node is taken by one vehicle
        model.constraint(model.partition(routes))

        quantities = model.array(quantities_data)
        time_warehouse = model.array(time_warehouse_data)
        time_matrix = model.array(time_matrix_data)
        loading_times = model.array(loading_times_data)
        starts = model.array(starts_data)
        ends = model.array(ends_data)
        # waiting[k] is the waiting time at node k
        waiting = [model.float(0, depot_tw_end) for k in range(nb_nodes)]
        waiting_array = model.array(waiting)
        distance_matrix = model.array(distance_matrix_data)
        distance_warehouse = model.array(distance_warehouse_data)

        times = [None] * nb_vehicles
        lateness = [None] * nb_vehicles
        home_lateness = [None] * nb_vehicles
        route_distances = [None] * nb_vehicles

        for k in range(nb_vehicles):
            route = routes[k]
            c = model.count(route)

            demand_lambda = model.lambda_function(lambda i, prev: prev + quantities[route[i]])
            # route_quantities[k][i] indicates the number of clients in vehicle k
            # at its i-th taken node
            route_quantities = model.array(model.range(0, c), demand_lambda)
            quantity_lambda = model.lambda_function(lambda i: route_quantities[i] <= capacity)
            # Vehicles have a maximum capacity
            model.constraint(model.and_(model.range(0, c), quantity_lambda))

            times_lambda = model.lambda_function(
                lambda i, prev: model.max(
                    starts[route[i]],
                    model.iif(
                        i == 0,
                        depot_starts[k] + time_warehouse[route[0]],
                        prev + time_matrix[route[i-1]][route[i]]
                    )
                ) + waiting_array[route[i]] + loading_times[route[i]]
            )
            # times[k][i] is the time at which vehicle k leaves the i-th node
            # (after waiting and loading time at node i)
            times[k] = model.array(model.range(0, c), times_lambda)

            lateness_lambda = model.lambda_function(
                lambda i: model.max(
                    0,
                    times[k][i] - loading_times[route[i]] - ends[route[i]]
                )
            )
            # Total lateness of the k-th route
            lateness[k] = model.sum(model.range(0, c), lateness_lambda)

            home_lateness[k] = model.iif(
                c > 0,
                model.max(0, times[k][c-1] + time_warehouse[route[c-1]] - depot_tw_end),
                0
            )

            route_dist_lambda = model.lambda_function(
                lambda i: distance_matrix[route[i-1]][route[i]]
            )
            route_distances[k] = model.sum(
                model.range(1, c),
                route_dist_lambda
            ) + model.iif(
                c > 0,
                distance_warehouse[route[0]] + distance_warehouse[route[c-1]],
                0
            )

        routes_array = model.array(routes)
        times_array = model.array(times)
        client_lateness = [None] * nb_clients
        for k in range(nb_clients):
            # For each pickup node k, its associated delivery node is k + nb_clients
            pickup_list_index = model.find(routes_array, k)
            delivery_list_index = model.find(routes_array, k + nb_clients)
            # A client picked up in route i is delivered in route i
            model.constraint(pickup_list_index == delivery_list_index)

            client_list = routes_array[pickup_list_index]
            pickup_index = model.index(client_list, k)
            delivery_index = model.index(client_list, k + nb_clients)
            # Pickup before delivery
            model.constraint(pickup_index < delivery_index)

            pickup_time = times_array[pickup_list_index][pickup_index]
            delivery_time = times_array[delivery_list_index][delivery_index] \
                - loading_times[k + nb_clients]
            travel_time = delivery_time - pickup_time
            client_lateness[k] = model.max(travel_time - max_travel_times[k], 0)

        total_lateness = model.sum(lateness + home_lateness)
        total_client_lateness = model.sum(client_lateness)
        total_distance = model.sum(route_distances)

        model.minimize(total_lateness)
        model.minimize(total_client_lateness)
        model.minimize(total_distance / scale)

        model.close()
        optimizer.param.time_limit = int(str_time_limit)
        optimizer.solve()

        #
        # Write the solution in a file with the following format:
        #  - total lateness on the routes, total client lateness, and total distance
        #  - for each vehicle, the depot start time, the nodes visited (omitting the start/end at the
        # depot), and the waiting time at each node
        #
        if sol_file is not None:
            with open(sol_file, 'w') as f:
                f.write("%d %d %.2f\n" % (
                    total_lateness.value,
                    total_client_lateness.value,
                    total_distance.value
                ))
                for k in range(nb_vehicles):
                    f.write("Vehicle %d (%.2f): " %(k + 1, depot_starts[k].value))
                    for node in routes[k].value:
                        f.write("%d (%.2f), " % (node, waiting[node].value))
                    f.write("\n")
    return 0

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

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

    main(instance_file, str_time_limit, sol_file)

#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <json.hpp>
#include <sstream>
#include <vector>

using namespace hexaly;
using namespace std;
using json = nlohmann::json;

class Darp {
public:
    HexalyOptimizer optimizer;

    int nbClients;
    int nbNodes;
    int nbVehicles;
    double factor;
    double depotTwEnd;
    int capacity;
    double scale;

    vector<int> quantitiesData;
    vector<vector<double>> distances;
    vector<double> startsData;
    vector<double> endsData;
    vector<double> loadingTimesData;
    vector<double> maxTravelTimes;
    vector<double> distanceWarehouseData;
    vector<double> timeWarehouseData;
    vector<vector<double>> distanceMatrixData;
    vector<vector<double>> timeMatrixData;

    vector<HxExpression> routes;
    vector<HxExpression> depotStarts;
    vector<HxExpression> waiting;

    HxExpression totalLateness;
    HxExpression totalClientLateness;
    HxExpression totalDistance;

    void readInstance(const string& fileName) {
        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        json instance;
        instance << infile;
        infile.close();

        nbClients = instance["nbClients"];
        nbNodes = instance["nbNodes"];
        nbVehicles = instance["nbVehicles"];
        depotTwEnd = instance["depot"]["twEnd"];
        capacity = instance["capacity"];
        scale = instance["scale"];

        quantitiesData.resize(2 * nbClients);
        startsData.resize(2 * nbClients);
        endsData.resize(2 * nbClients);
        loadingTimesData.resize(2 * nbClients);
        maxTravelTimes.resize(2 * nbClients);

        for (int k = 0; k < nbClients; ++k) {
            quantitiesData[k] = instance["clients"][k]["nbClients"].get<int>();
            quantitiesData[k + nbClients] = -instance["clients"][k]["nbClients"].get<int>();

            startsData[k] = instance["clients"][k]["pickup"]["start"].get<double>();
            startsData[k + nbClients] = instance["clients"][k]["delivery"]["start"].get<double>();

            endsData[k] = instance["clients"][k]["pickup"]["end"].get<double>();
            endsData[k + nbClients] = instance["clients"][k]["delivery"]["end"].get<double>();

            loadingTimesData[k] = instance["clients"][k]["pickup"]["loadingTime"].get<double>();
            loadingTimesData[k + nbClients] = instance["clients"][k]["delivery"]["loadingTime"].get<double>();

            maxTravelTimes[k] = instance["clients"][k]["pickup"]["maxTravelTime"].get<double>();
            maxTravelTimes[k + nbClients] = instance["clients"][k]["delivery"]["maxTravelTime"].get<double>();
        }

        distances.resize(nbNodes + 1, vector<double>(nbNodes + 1));
        for (int i = 0; i < nbNodes + 1; ++i) {
            for (int j = 0; j < nbNodes + 1; ++j) {
                distances[i][j] = instance["distanceMatrix"][i][j].get<double>();
            }
        }

        factor = 1.0 / (instance["scale"].get<double>() * instance["speed"].get<double>());

        distanceWarehouseData.resize(nbNodes);
        timeWarehouseData.resize(nbNodes);
        for (int k = 0; k < nbNodes; ++k) {
            distanceWarehouseData[k] = distances[0][k+1];
            timeWarehouseData[k] = distanceWarehouseData[k] * factor;
        }

        distanceMatrixData.resize(nbNodes, vector<double>(nbNodes));
        timeMatrixData.resize(nbNodes, vector<double>(nbNodes));
        for (int i = 0; i < nbNodes; ++i) {
            for (int j = 0; j < nbNodes; ++j){
                distanceMatrixData[i][j] = distances[i+1][j+1];
                timeMatrixData[i][j] = distanceMatrixData[i][j] * factor;
            }
        }
    }
    
    void solve(int limit) {
        HxModel model = optimizer.getModel();

        // routes[k] represents the nodes visited by vehicle k
        routes.resize(nbVehicles);
        depotStarts.resize(nbVehicles);
        for (int k = 0; k < nbVehicles; ++k) {
            routes[k] = model.listVar(nbNodes);
            depotStarts[k] = model.floatVar(0.0, depotTwEnd);
        }
        // waiting[k] is the waiting time at node k
        waiting.resize(nbNodes);
        for (int k = 0; k < nbNodes; ++k) {
            waiting[k] = model.floatVar(0.0, depotTwEnd);
        }
        // Each node is taken by one vehicle
        model.constraint(model.partition(routes.begin(), routes.end()));

        HxExpression quantities = model.array(quantitiesData.begin(), quantitiesData.end());
        HxExpression timeWarehouse = model.array(timeWarehouseData.begin(), timeWarehouseData.end());
        HxExpression timeMatrix = model.array();
        for (int i = 0; i < nbNodes; ++i) {
            timeMatrix.addOperand(model.array(timeMatrixData[i].begin(), timeMatrixData[i].end()));
        }
        HxExpression loadingTimes = model.array(loadingTimesData.begin(), loadingTimesData.end());
        HxExpression starts = model.array(startsData.begin(), startsData.end());
        HxExpression ends = model.array(endsData.begin(), endsData.end());
        HxExpression waitingArray = model.array(waiting.begin(), waiting.end());
        HxExpression distanceMatrix = model.array();
        for (int i = 0; i < nbNodes; ++i) {
            distanceMatrix.addOperand(model.array(distanceMatrixData[i].begin(), distanceMatrixData[i].end()));
        }
        HxExpression distanceWarehouse = model.array(distanceWarehouseData.begin(), distanceWarehouseData.end());
        vector<HxExpression> times(nbVehicles);
        vector<HxExpression> lateness(nbVehicles);
        vector<HxExpression> homeLateness(nbVehicles);
        vector<HxExpression> routeDistances(nbVehicles);

        for (int k = 0; k < nbVehicles; ++k) {
            HxExpression route = routes[k];
            HxExpression c = model.count(route);

            HxExpression demandLambda = model.createLambdaFunction(
                [&](HxExpression i, HxExpression prev) { return prev + quantities[route[i]]; });
            // routeQuantities[k][i] indicates the number of clients in vehicle k
            // at its i-th taken node
            HxExpression routeQuantities = model.array(model.range(0, c), demandLambda);
            HxExpression quantityLambda = model.createLambdaFunction(
                [&](HxExpression i) { return (routeQuantities[i] <= capacity); });
            // Vehicles have a maximum capacity
            model.constraint(model.and_(model.range(0, c), quantityLambda));

            HxExpression timesLambda = model.createLambdaFunction(
                [&](HxExpression i, HxExpression prev) {
                    return model.max(starts[route[i]], model.iif(
                        i == 0,
                        depotStarts[k] + timeWarehouse[route[0]],
                        prev + timeMatrix[route[i-1]][route[i]]
                    )) + waitingArray[route[i]] + loadingTimes[route[i]];
                }
            );
            // times[k][i] is the time at which vehicle k leaves the i-th node
            // (after waiting and loading time at node i)
            times[k] = model.array(model.range(0, c), timesLambda);

            HxExpression latenessLambda = model.createLambdaFunction(
                [&](HxExpression i) {
                    return model.max(
                        0,
                        times[k][i] - loadingTimes[route[i]] - ends[route[i]]
                    );
                }
            );
            // Total lateness of the k-th route
            lateness[k] = model.sum(model.range(0, c), latenessLambda);

            homeLateness[k] = model.iif(
                c > 0,
                model.max(0, times[k][c-1] + timeWarehouse[route[c-1]] - depotTwEnd),
                0
            );

            HxExpression routeDistLambda = model.createLambdaFunction(
                [&](HxExpression i) { return distanceMatrix[route[i-1]][route[i]]; });
            routeDistances[k] = model.sum(model.range(1, c), routeDistLambda)
                + model.iif(
                    c > 0,
                    distanceWarehouse[route[0]] + distanceWarehouse[route[c-1]],
                    0
                );
        }

        HxExpression routesArray = model.array(routes.begin(), routes.end());
        HxExpression timesArray = model.array(times.begin(), times.end());
        vector<HxExpression> clientLateness(nbClients);

        for (int k = 0; k < nbClients; ++k) {
            // For each pickup node k, its associated delivery node is k + instance.nbClients
            HxExpression pickupListIndex = model.find(routesArray, k);
            HxExpression deliveryListIndex = model.find(routesArray, k + nbClients);
            // A client picked up in route i is delivered in route i
            model.constraint(pickupListIndex == deliveryListIndex);

            HxExpression clientList = routesArray[pickupListIndex];
            HxExpression pickupIndex = model.indexOf(clientList, k);
            HxExpression deliveryIndex = model.indexOf(clientList, k + nbClients);
            // Pickup before delivery
            model.constraint(pickupIndex < deliveryIndex);

            HxExpression pickupTime = timesArray[pickupListIndex][pickupIndex];
            HxExpression deliveryTime =
                timesArray[deliveryListIndex][deliveryIndex] - loadingTimes[k + nbClients];
            HxExpression travelTime = deliveryTime - pickupTime;
            clientLateness[k] = model.max(travelTime - maxTravelTimes[k], 0);
        }

        vector<HxExpression> latenessPlusHomeLateness(nbVehicles);
        for (int k = 0; k < nbVehicles; ++k) {
            latenessPlusHomeLateness[k] = lateness[k] + homeLateness[k];
        }

        totalLateness = model.sum(latenessPlusHomeLateness.begin(), latenessPlusHomeLateness.end());
        totalClientLateness = model.sum(clientLateness.begin(), clientLateness.end());
        totalDistance = model.sum(routeDistances.begin(), routeDistances.end());

        model.minimize(totalLateness);
        model.minimize(totalClientLateness);
        model.minimize(totalDistance / scale);

        model.close();

        optimizer.getParam().setTimeLimit(limit);

        optimizer.getParam().setSeed(5);
        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - total lateness on the routes, total client lateness, total distance
     *  - for each vehicle, the depot start time, the nodes visited (omitting the start/end at the
     * depot), and the waiting time at each node */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << totalLateness.getDoubleValue() << " " << totalClientLateness.getDoubleValue()
            << " " << totalDistance.getDoubleValue() << endl;
        for (int k = 0; k < nbVehicles; ++k) {
            HxCollection route = routes[k].getCollectionValue();
            outfile << "Vehicle " << k << " (" << depotStarts[k].getDoubleValue() << "): ";
            for (int i = 0; i < route.count(); ++i) {
                outfile << route[i] << " (" << waiting[route[i]].getDoubleValue() << "), ";
            }
            outfile << endl;
        }
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: darp 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 {
        Darp 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;
using Newtonsoft.Json.Linq;
using Newtonsoft.Json;

public class Darp : IDisposable
{
    HexalyOptimizer optimizer;

    int nbClients;
    int nbNodes;
    int nbVehicles;
    double depotTwEnd;
    int capacity;
    double scale;
    double factor;

    int[] quantitiesData;
    double[] startsData;
    double[] endsData;
    double[] loadingTimesData;
    double[] maxTravelTimes;

    double[][] distances;
    double[] distanceWarehouseData;
    double[] timeWarehouseData;
    double[][] distanceMatrixData;
    double[][] timeMatrixData;

    HxExpression[] routes;
    HxExpression[] depotStarts;
    HxExpression[] waiting;

    HxExpression[] clientLateness;
    HxExpression[] latenessPlusHomeLateness;

    HxExpression totalLateness;
    HxExpression totalClientLateness;
    HxExpression totalDistance;

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

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        JObject instance = JObject.Parse(File.ReadAllText(fileName));

        nbClients = (int)instance["nbClients"];
        nbNodes = (int)instance["nbNodes"];
        nbVehicles = (int)instance["nbVehicles"];
        depotTwEnd = (double)instance["depot"]["twEnd"];
        capacity = (int)instance["capacity"];
        scale = (double)instance["scale"];

        quantitiesData = new int[2 * nbClients];
        startsData = new double[2 * nbClients];
        endsData = new double[2 * nbClients];
        loadingTimesData = new double[2 * nbClients];
        maxTravelTimes = new double[2 * nbClients];

        for (int k = 0; k < nbClients; ++k)
        {
            quantitiesData[k] = (int)instance["clients"][k]["nbClients"];
            quantitiesData[k + nbClients] = -(int)instance["clients"][k]["nbClients"];

            startsData[k] = (double)instance["clients"][k]["pickup"]["start"];
            startsData[k + nbClients] = (double)instance["clients"][k]["delivery"]["start"];

            endsData[k] = (double)instance["clients"][k]["pickup"]["end"];
            endsData[k + nbClients] = (double)instance["clients"][k]["delivery"]["end"];

            loadingTimesData[k] = (double)instance["clients"][k]["pickup"]["loadingTime"];
            loadingTimesData[k + nbClients] = (double)instance["clients"][k]["delivery"]["loadingTime"];

            maxTravelTimes[k] = (double)instance["clients"][k]["pickup"]["maxTravelTime"];
            maxTravelTimes[k + nbClients] =
                (double)instance["clients"][k]["delivery"]["maxTravelTime"];
        }

        distances = new double[nbNodes + 1][];
        for (int i = 0; i < nbNodes + 1; ++i)
        {
            distances[i] = new double[nbNodes + 1];
            for (int j = 0; j < nbNodes + 1; ++j)
            {
                distances[i][j] = (double)instance["distanceMatrix"][i][j];
            }
        }

        factor = 1 / ((double)instance["scale"] * (double)instance["speed"]);

        distanceWarehouseData = new double[nbNodes];
        timeWarehouseData = new double[nbNodes];
        for (int k = 0; k < nbNodes; ++k)
        {
            distanceWarehouseData[k] = distances[0][k+1];
            timeWarehouseData[k] = distanceWarehouseData[k] * factor;
        }

        distanceMatrixData = new double[nbNodes][];
        timeMatrixData = new double[nbNodes][];
        for (int i = 0; i < nbNodes; ++i)
        {
            distanceMatrixData[i] = new double[nbNodes];
            timeMatrixData[i] = new double[nbNodes];
            for (int j = 0; j < nbNodes; ++j)
            {
                distanceMatrixData[i][j] = distances[i+1][j+1];
                timeMatrixData[i][j] = distanceMatrixData[i][j] * factor;
            }
        }

    }

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

    void Solve(int limit)
    {
        HxModel model = optimizer.GetModel();

        // routes[k] represents the nodes visited by vehicle k
        routes = new HxExpression[nbVehicles];
        depotStarts = new HxExpression[nbVehicles];
        for (int k = 0; k < nbVehicles; ++k)
        {
            routes[k] = model.List(nbNodes);
            depotStarts[k] = model.Float(0.0, depotTwEnd);
        }
        // waiting[k] is the waiting time at node k
        waiting = new HxExpression[nbNodes];
        for (int k = 0; k < nbNodes; ++k)
        {
            waiting[k] = model.Float(0.0, depotTwEnd);
        }
        // Each node is taken by one vehicle
        model.Constraint(model.Partition(routes));

        HxExpression quantities = model.Array(quantitiesData);
        HxExpression timeWarehouse = model.Array(timeWarehouseData);
        HxExpression loadingTimes = model.Array(loadingTimesData);
        HxExpression starts = model.Array(startsData);
        HxExpression ends = model.Array(endsData);
        HxExpression waitingArray = model.Array(waiting);
        HxExpression timeMatrix = model.Array();
        HxExpression distanceMatrix = model.Array();
        for (int i = 0; i < nbNodes; ++i)
        {
            timeMatrix.AddOperand(model.Array(timeMatrixData[i]));
            distanceMatrix.AddOperand(model.Array(distanceMatrixData[i]));
        }
        HxExpression distanceWarehouse = model.Array(distanceWarehouseData);

        HxExpression[] times = new HxExpression[nbVehicles];
        HxExpression[] lateness = new HxExpression[nbVehicles];
        HxExpression[] homeLateness = new HxExpression[nbVehicles];
        HxExpression[] routeDistances = new HxExpression[nbVehicles];

        for (int k = 0; k < nbVehicles; ++k)
        {
            HxExpression route = routes[k];
            HxExpression c = model.Count(route);

            HxExpression demandLambda = model.LambdaFunction(
                (i, prev) => prev + quantities[route[i]]
            );
            // routeQuantities[k][i] indicates the number of clients in vehicle k
            // at its i-th taken node
            HxExpression routeQuantities = model.Array(model.Range(0, c), demandLambda);
            HxExpression quantityLambda = model.LambdaFunction(i => routeQuantities[i] <= capacity);
            // Vehicles have a maximum capacity
            model.Constraint(model.And(model.Range(0, c), quantityLambda));

            HxExpression timesLambda = model.LambdaFunction(
                (i, prev) =>
                    model.Max(
                        starts[route[i]],
                        model.If(
                            i == 0,
                            depotStarts[k] + timeWarehouse[route[0]],
                            prev + timeMatrix[route[i-1]][route[i]]
                        )
                    ) + waitingArray[route[i]] + loadingTimes[route[i]]);
            // times[k][i] is the time at which vehicle k leaves the i-th node
            // (after waiting and loading time at node i)
            times[k] = model.Array(model.Range(0, c), timesLambda);

            HxExpression latenessLambda = model.LambdaFunction(
                i =>
                    model.Max(
                        0,
                        times[k][i] - loadingTimes[route[i]] - ends[route[i]]
                    )
            );
            // Total lateness of the k-th route
            lateness[k] = model.Sum(model.Range(0, c), latenessLambda);

            homeLateness[k] = model.If(
                c > 0,
                model.Max(0, times[k][c-1] + timeWarehouse[route[c-1]] - depotTwEnd),
                0
            );

            HxExpression routeDistLambda = model.LambdaFunction(
                i => distanceMatrix[route[i-1]][route[i]]
            );
            routeDistances[k] = model.Sum(
                model.Range(1, c),
                routeDistLambda
            ) + model.If(
                c > 0,
                distanceWarehouse[route[0]] + distanceWarehouse[route[c-1]],
                0
            );
        }
     
        HxExpression routesArray = model.Array(routes);
        HxExpression timesArray = model.Array(times);
        clientLateness = new HxExpression[nbClients];

        for (int k = 0; k < nbClients; ++k)
        {
            // For each pickup node k, its associated delivery node is k + nbClients
            HxExpression pickupListIndex = model.Find(routesArray, k);
            HxExpression deliveryListIndex = model.Find(routesArray, k + nbClients);
            // A client picked up in route i is delivered in route i
            model.Constraint(pickupListIndex == deliveryListIndex);

            HxExpression clientList = routesArray[pickupListIndex];
            HxExpression pickupIndex = model.IndexOf(clientList, k);
            HxExpression deliveryIndex = model.IndexOf(clientList, k + nbClients);
            // Pickup before delivery
            model.Constraint(pickupIndex < deliveryIndex);

            HxExpression pickupTime = timesArray[pickupListIndex][pickupIndex];
            HxExpression deliveryTime =
                timesArray[deliveryListIndex][deliveryIndex] - loadingTimes[k + nbClients];
            HxExpression travelTime = deliveryTime - pickupTime;
            clientLateness[k] = model.Max(travelTime - maxTravelTimes[k], 0);
        }

        latenessPlusHomeLateness = new HxExpression[nbVehicles];
        for (int k = 0; k < nbVehicles; ++k)
        {
            latenessPlusHomeLateness[k] = lateness[k] + homeLateness[k];
        }
        totalLateness = model.Sum(latenessPlusHomeLateness);
        totalClientLateness = model.Sum(clientLateness);
        totalDistance = model.Sum(routeDistances);

        model.Minimize(totalLateness);
        model.Minimize(totalClientLateness);
        model.Minimize(totalDistance / scale);

        model.Close();

        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.GetParam().SetSeed(5);
        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - total lateness on the routes, total client lateness, total distance
     *  - for each vehicle, the depot start time, the nodes visited (omitting the start/end at the
     * depot), and the waiting time at each node */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(
                totalLateness.GetDoubleValue()
                + " "
                + totalClientLateness.GetDoubleValue()
                + " "
                + totalDistance.GetDoubleValue()
            );
            for (int k = 0; k < nbVehicles; ++k)
            {
                HxCollection route = routes[k].GetCollectionValue();
                output.Write("Vehicle " + (k + 1) + " (" + depotStarts[k].GetDoubleValue() + "): ");
                for (int i = 0; i < route.Count(); ++i)
                    output.Write(route[i] + " (" + (waiting[route[i]].GetDoubleValue() + "), "));
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Darp 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 (Darp model = new Darp())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}

import com.hexaly.optimizer.*;

import java.io.*;
import java.util.*;

import com.google.gson.JsonObject;
import com.google.gson.JsonArray;
import com.google.gson.JsonParser;
import com.google.gson.stream.JsonReader;

public class Darp {
    private final HexalyOptimizer optimizer;

    int nbClients;
    int nbNodes;
    int nbVehicles;
    double depotTwEnd;
    int capacity;
    double scale;
    double factor;

    int[] quantitiesData;
    double[] startsData;
    double[] endsData;
    double[] loadingTimesData;
    double[] maxTravelTimes;

    double[][] distances;
    double[] distanceWarehouseData;
    double[] timeWarehouseData;
    double[][] distanceMatrixData;
    double[][] timeMatrixData;

    HxExpression[] routes;
    HxExpression[] depotStarts;
    HxExpression[] waiting;

    HxExpression totalLateness;
    HxExpression totalClientLateness;
    HxExpression totalDistance;

    private Darp() {
        optimizer = new HexalyOptimizer();
    }

    // Read instance data
    private void readInstance(String fileName) throws IOException {
        JsonReader reader = new JsonReader(new InputStreamReader(new FileInputStream(fileName)));
        JsonObject instance = JsonParser.parseReader(reader).getAsJsonObject();

        nbClients = instance.get("nbClients").getAsInt();
        nbNodes = instance.get("nbNodes").getAsInt();
        nbVehicles = instance.get("nbVehicles").getAsInt();
        JsonObject depot = (JsonObject) instance.get("depot");
        depotTwEnd = depot.get("twEnd").getAsDouble();
        capacity = instance.get("capacity").getAsInt();
        scale = instance.get("scale").getAsDouble();
        JsonArray clients = (JsonArray) instance.get("clients");

        quantitiesData = new int[2 * nbClients];
        startsData = new double[2 * nbClients];
        endsData = new double[2 * nbClients];
        loadingTimesData = new double[2 * nbClients];
        maxTravelTimes = new double[2 * nbClients];

        for (int k = 0; k < nbClients; ++k){
            JsonObject clientk = (JsonObject) clients.get(k);

            quantitiesData[k] = clientk.get("nbClients").getAsInt();
            quantitiesData[k + nbClients] = -clientk.get("nbClients").getAsInt();

            JsonObject clientkPickup = (JsonObject) clientk.get("pickup");
            JsonObject clientkDelivery = (JsonObject) clientk.get("delivery");

            startsData[k] = clientkPickup.get("start").getAsDouble();
            startsData[k + nbClients] = clientkDelivery.get("start").getAsDouble();

            endsData[k] = clientkPickup.get("end").getAsDouble();
            endsData[k + nbClients] = clientkDelivery.get("end").getAsDouble();

            loadingTimesData[k] = clientkPickup.get("loadingTime").getAsDouble();
            loadingTimesData[k + nbClients] = clientkDelivery.get("loadingTime").getAsDouble();

            maxTravelTimes[k] = clientkPickup.get("maxTravelTime").getAsDouble();
            maxTravelTimes[k + nbClients] = clientkDelivery.get("maxTravelTime").getAsDouble();
        }

        distances = new double[nbNodes + 1][nbNodes + 1];
        JsonArray distanceMatrixJson = (JsonArray) instance.get("distanceMatrix");
        for (int i = 0; i < nbNodes + 1; ++i){
            JsonArray distanceMatrixJsoni = (JsonArray) distanceMatrixJson.get(i);
            for (int j = 0; j < nbNodes + 1; ++j){
                distances[i][j] = distanceMatrixJsoni.get(j).getAsDouble();
            }
        }

        factor = 1.0 / (instance.get("scale").getAsDouble() * instance.get("speed").getAsDouble());

        distanceWarehouseData = new double[nbNodes];
        timeWarehouseData = new double[nbNodes];
        for (int k = 0; k < nbNodes; ++k){
            distanceWarehouseData[k] = distances[0][k+1];
            timeWarehouseData[k] = distanceWarehouseData[k] * factor;
        }

        distanceMatrixData = new double[nbNodes][nbNodes];
        timeMatrixData = new double[nbNodes][nbNodes];
        for (int i = 0; i < nbNodes; ++i){
            for (int j = 0; j < nbNodes; ++j){
                distanceMatrixData[i][j] = distances[i+1][j+1];
                timeMatrixData[i][j] = distanceMatrixData[i][j] * factor;
            }
        }
    }

    private void solve(int limit) {
        HxModel model = optimizer.getModel();

        // routes[k] represents the nodes visited by vehicle k
        routes = new HxExpression[nbVehicles];
        depotStarts = new HxExpression[nbVehicles];
        for (int k = 0; k < nbVehicles; ++k){
            routes[k] = model.listVar(nbNodes);
            depotStarts[k] = model.floatVar(0.0, depotTwEnd);
        }
        // waiting[k] is the waiting time at node k
        waiting = new HxExpression[nbNodes];
        for (int k = 0; k < nbNodes; ++k){
            waiting[k] = model.floatVar(0.0, depotTwEnd);
        }
        // Each node is taken by one vehicle
        model.constraint(model.partition(routes));

        HxExpression quantities = model.array(quantitiesData);
        HxExpression timeWarehouse = model.array(timeWarehouseData);
        HxExpression timeMatrix = model.array(timeMatrixData);
        for (int i = 0; i < nbNodes; ++i){
            timeMatrix.addOperand(model.array(timeMatrixData[i]));
        }
        HxExpression loadingTimes = model.array(loadingTimesData);
        HxExpression starts = model.array(startsData);
        HxExpression ends = model.array(endsData);
        HxExpression waitingArray = model.array(waiting);
        HxExpression distanceMatrix = model.array();
        for (int i = 0; i < nbNodes; ++i){
            distanceMatrix.addOperand(model.array(distanceMatrixData[i]));
        }
        HxExpression distanceWarehouse = model.array(distanceWarehouseData);

        HxExpression[] times = new HxExpression[nbVehicles];
        HxExpression[] lateness = new HxExpression[nbVehicles];
        HxExpression[] homeLateness = new HxExpression[nbVehicles];
        HxExpression[] routeDistances = new HxExpression[nbVehicles];

        for (int k = 0; k < nbVehicles; ++k){
            HxExpression route = routes[k];
            HxExpression c = model.count(route);

            HxExpression demandLambda = model.lambdaFunction(
                (i, prev) -> model.sum(prev, model.at(quantities, model.at(route, i))));
            // routeQuantities[k][i] indicates the number of clients in vehicle k
            // at its i-th taken node
            HxExpression routeQuantities = model.array(model.range(0, c), demandLambda);
            HxExpression quantityLambda = model.createLambdaFunction(
                i -> model.leq(model.at(routeQuantities, i), capacity));
            // Vehicles have a maximum capacity
            model.constraint(model.and(model.range(0, c), quantityLambda));

            HxExpression depotStartsk = depotStarts[k];
            HxExpression timesLambda = model.lambdaFunction(
                (i, prev) -> model.sum(
                    model.max(
                        model.at(starts, model.at(route, i)),
                        model.iif(
                            model.eq(i, 0),
                            model.sum(
                                depotStartsk,
                                model.at(timeWarehouse, model.at(route, 0))
                            ),
                            model.sum(
                                prev,
                                model.at(
                                    timeMatrix,
                                    model.at(route, model.sub(i, 1)),
                                    model.at(route, i)
                                )
                            )
                        )
                    ),
                    model.sum(
                        model.at(waitingArray, model.at(route, i)),
                        model.at(loadingTimes, model.at(route, i))
                    )
                )
            );
            HxExpression timesk = model.array(model.range(0, c), timesLambda);
            // times[k][i] is the time at which vehicle k leaves the i-th node
            // (after waiting and loading time at node i)
            times[k] = timesk;
            HxExpression latenessLambda = model.lambdaFunction(
                i -> model.max(
                    0,
                    model.sub(
                        model.sub(
                            model.at(timesk, i),
                            model.at(loadingTimes, model.at(route, i))
                        ),
                        model.at(ends, model.at(route, i))
                    )
                )
            );
            // Total lateness of the k-th route
            lateness[k] = model.sum(model.range(0, c), latenessLambda);

            homeLateness[k] = model.iif(
                model.gt(c, 0),
                model.max(
                    0,
                    model.sum(
                        model.at(timesk, model.sub(c, 1)),
                        model.sub(
                            model.at(timeWarehouse, model.at(route, model.sub(c, 1))),
                            depotTwEnd
                        )
                    )
                ),
                0
            );

            HxExpression routeDistLambda = model.lambdaFunction(
                i -> model.at(
                    distanceMatrix,
                    model.at(route, model.sub(i, 1)),
                    model.at(route, i)
                )
            );
            routeDistances[k] = model.sum(
                model.sum(model.range(1, c), routeDistLambda),
                model.iif(
                    model.gt(c, 0),
                    model.sum(
                        model.at(distanceWarehouse, model.at(route, 0)),
                        model.at(distanceWarehouse, model.at(route, model.sub(c, 1)))
                    ),
                    0
                )
            );
        }

        HxExpression routesArray = model.array(routes);
        HxExpression timesArray = model.array(times);
        HxExpression[] clientLateness = new HxExpression[nbClients];

        for (int k = 0; k < nbClients; ++k) {
            // For each pickup node k, its associated delivery node is k + nbClients
            HxExpression pickupListIndex = model.find(routesArray, k);
            HxExpression deliveryListIndex = model.find(routesArray, k + nbClients);
            // A client picked up in route i is delivered in route i
            model.constraint(model.eq(pickupListIndex, deliveryListIndex));

            HxExpression clientList = model.at(routesArray, pickupListIndex);
            HxExpression pickupIndex = model.indexOf(clientList, k);
            HxExpression deliveryIndex = model.indexOf(clientList, k + nbClients);
            // Pickup before delivery
            model.constraint(model.lt(pickupIndex, deliveryIndex));

            HxExpression pickupTime = model.at(timesArray, pickupListIndex, pickupIndex);
            HxExpression deliveryTime = model.sub(
                model.at(timesArray, deliveryListIndex, deliveryIndex),
                model.at(loadingTimes, k + nbClients)
            );
            HxExpression travelTime = model.sub(deliveryTime, pickupTime);
            clientLateness[k] = model.max(model.sub(travelTime, maxTravelTimes[k]), 0);
        }

        HxExpression[] latenessPlusHomeLateness = new HxExpression[nbVehicles];
        for (int k = 0; k < nbVehicles; ++k) {
            latenessPlusHomeLateness[k] = model.sum(lateness[k], homeLateness[k]);
        }

        totalLateness = model.sum(latenessPlusHomeLateness);
        totalClientLateness = model.sum(clientLateness);
        totalDistance = model.sum(routeDistances);

        model.minimize(totalLateness);
        model.minimize(totalClientLateness);
        model.minimize(model.div(totalDistance, scale));

        model.close();

        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - total lateness on the routes, total client lateness, total distance
     *  - for each vehicle, the depot start time, the nodes visited (omitting the start/end at the
     * depot), and the waiting time at each node */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(
                totalLateness.getDoubleValue()
                + " "
                + totalClientLateness.getDoubleValue()
                + " "
                + totalDistance.getDoubleValue()
            );
            for (int k = 0; k < nbVehicles; ++k) {
                HxCollection route = routes[k].getCollectionValue();
                output.print("Vehicle " + (k + 1) + " (" + depotStarts[k].getDoubleValue() + "): ");
                for (int i = 0; i < route.count(); ++i) {
                    int routei = (int) route.get(i); // route.get(i) casted to long
                    output.print(routei + " (" + (waiting[routei].getDoubleValue() + "), "));
                }
                output.println();
            }
        }
    }

    public static void main(String[] args) {

        if (args.length < 1) {
            System.err.println("Usage: java Darp inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

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

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