Split Delivery Vehicle Routing Problem (SDVRP)

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, sol_file):
    #
    # Read instance data
    #
    nb_customers, truck_capacity, dist_matrix_data, dist_depot_data, demands_data = \
        read_input_sdvrp(instance_file)
    nb_trucks = nb_customers

    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 _ in range(nb_trucks)]

        # Quantity carried by each truck for each customer
        quantity = [None] * nb_trucks
        for k in range(nb_trucks):
            quantity[k] = [model.float(0, demands_data[i]) for i in range(nb_customers)]

        # All customers must be visited at least by one truck
        model.constraint(model.cover(customers_sequences))

        # Create Hexaly arrays to be able to access them with "at" operators
        dist_matrix = model.array(dist_matrix_data)
        dist_depot = model.array(dist_depot_data)

        route_distances = [None] * nb_trucks
        trucks_used = [None] * nb_trucks

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

            # A truck is used if it visits at least one customer
            trucks_used[k] = model.count(sequence) > 0

            # The quantity carried in each route must not exceed the truck capacity
            quantity_array = model.array(quantity[k])
            quantity_lambda = model.lambda_function(lambda j: quantity_array[j])
            route_quantity = model.sum(sequence, quantity_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]))
            route_distances[k] = model.sum(model.range(1, c), dist_lambda) \
                + model.iif(
                    trucks_used[k],
                    dist_depot[sequence[0]] + dist_depot[sequence[c - 1]],
                    0)

        for i in range(nb_customers):
            # Each customer must receive at least its demand
            quantity_served = model.sum(
                quantity[k][i] * model.contains(customers_sequences[k], i)
                for k in range(nb_trucks))
            model.constraint(quantity_served >= demands_data[i])

        total_distance = model.sum(route_distances)

        # Objective: minimize the distance traveled
        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:
        # each line k contain the customers visited by the truck k
        #
        if len(sys.argv) >= 3:
            with open(sol_file, 'w') as f:
                f.write("%d \n" % 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
                    for customer in customers_sequences[k].value:
                        f.write("%d " % (customer + 1))
                    f.write("\n")


def read_input_sdvrp(filename):
    file_it = iter(read_elem(filename))

    nb_customers = int(next(file_it))
    capacity = int(next(file_it))

    demands = [None] * nb_customers
    for i in range(nb_customers):
        demands[i] = int(next(file_it))

    # Extracting the coordinates of the depot and the customers
    customers_x = [None] * nb_customers
    customers_y = [None] * nb_customers
    depot_x = float(next(file_it))
    depot_y = float(next(file_it))
    for i in range(nb_customers):
        customers_x[i] = float(next(file_it))
        customers_y[i] = float(next(file_it))

    distance_matrix = compute_distance_matrix(customers_x, customers_y)
    distance_depots = compute_distance_depots(depot_x, depot_y, customers_x, customers_y)
    return nb_customers, capacity, distance_matrix, distance_depots, demands


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


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


def compute_dist(xi, xj, yi, yj):
    exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
    return int(math.floor(exact_dist + 0.5))


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python sdvrp.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 <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand of each customer
    vector<int> demandsData;

    // Distance matrix
    vector<vector<int>> distMatrixData;

    // Distance between customers and depot
    vector<int> distDepotData;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    vector<HxExpression> customersSequences;

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    // Quantity carried by each truck
    vector<vector<HxExpression>> quantity;

    /* Read instance data */
    void readInstance(const string& fileName) {
        readInputSdvrp(fileName);
        nbTrucks = nbCustomers;
    }

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

        // Sequence of customers visited by each truck
        customersSequences.resize(nbCustomers);
        for (int k = 0; k < nbTrucks; ++k) {
            customersSequences[k] = model.listVar(nbCustomers);
        }

        // Quantity carried by each truck for each customer
        quantity.resize(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            quantity[k].resize(nbTrucks);
            for (int i = 0; i < nbCustomers; ++i) {
                quantity[k][i] = model.floatVar(0, demandsData[i]);
            }
        }

        // All customers must be visited by exactly one truck
        model.constraint(model.cover(customersSequences.begin(), customersSequences.end()));

        // Create Hexaly arrays to be able to access them with "at" operators
        HxExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            distMatrix.addOperand(model.array(distMatrixData[n].begin(), distMatrixData[n].end()));
        }
        HxExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());

        trucksUsed.resize(nbTrucks);
        vector<HxExpression> distRoutes(nbTrucks);

        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 carried in each route must not exceed the truck capacity
            HxExpression quantityArray = model.array(quantity[k].begin(), quantity[k].end());
            HxExpression quantityLambda = model.createLambdaFunction(
                    [&](HxExpression j) { return quantityArray[j]; });
            HxExpression routeQuantity = model.sum(sequence, quantityLambda);
            model.constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            HxExpression distLambda = model.createLambdaFunction([&](HxExpression i) {
                return model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i));
            });
            distRoutes[k] = model.sum(model.sum(model.range(1, c), distLambda),
                                      model.iif(trucksUsed[k],
                                                model.sum(model.at(distDepot, model.at(sequence, 0)),
                                                          model.at(distDepot, model.at(sequence, c - 1))),
                                                0));
        }

        for (int i = 0; i < nbCustomers; ++i) {
            // Each customer must receive at least its demand
            HxExpression quantityServed = model.sum();
            for (int k = 0; k < nbTrucks; ++k)
                quantityServed.addOperand(quantity[k][i] * model.contains(customersSequences[k], i));
            model.constraint(quantityServed >= demandsData[i]);
        }

        totalDistance = model.sum(distRoutes.begin(), distRoutes.end());

        // Objective: minimize the distance traveled
        model.minimize(totalDistance);

        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - Each line k contains the customers visited by the truck k */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << totalDistance.getValue() << 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
            HxCollection customersCollection = customersSequences[k].getCollectionValue();
            for (int i = 0; i < customersCollection.count(); ++i) {
                outfile << customersCollection[i] + 1 << " ";
            }
            outfile << endl;
        }
    }

private:
    void readInputSdvrp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }

        string str;
        char* line;
        istringstream iss;
        getline(infile, str);
        line = strdup(str.c_str());
        iss.str(line);
        iss >> nbCustomers;
        iss >> truckCapacity;
        nbTrucks = nbCustomers;
        getline(infile, str);
        line = strdup(str.c_str());
        istringstream stream(line);
        demandsData.resize(nbCustomers);
        int demand;
        for (int n = 0; n < nbCustomers; ++n) {
            stream >> demand;
            demandsData[n] = demand;
        }
        vector<int> customersX(nbCustomers);
        vector<int> customersY(nbCustomers);
        int depotX, depotY;
        getline(infile, str);
        line = strdup(str.c_str());
        istringstream iss_1(line);
        iss_1 >> depotX;
        iss_1 >> depotY;
        for (int i = 0; i < nbCustomers; ++i) {
            getline(infile, str);
            line = strdup(str.c_str());
            istringstream iss(line);
            iss >> customersX[i];
            iss >> customersY[i];
        }

        computeDistanceMatrix(depotX, depotY, customersX, customersY);

        infile.close();
    }

    // Compute the distance matrix
    void computeDistanceMatrix(int depotX, int depotY, const vector<int>& customersX, const vector<int>& customersY) {
        distMatrixData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i].resize(nbCustomers);
        }
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                int distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = distance;
                distMatrixData[j][i] = distance;
            }
        }

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

    static int computeDist(int xi, int xj, int yi, int yj) {
        double exactDist = sqrt(pow((double)xi - xj, 2) + pow((double)yi - yj, 2));
        return floor(exactDist + 0.5);
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: Sdvrp 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 {
        Sdvrp 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;
    }
}

/********** Sdvrp.cs **********/
using System;
using System.IO;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand of each customer
    long[] demandsData;

    // Distance matrix
    long[][] distMatrixData;

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

    // Number of trucks
    int nbTrucks;

    // Decision variables
    HxExpression[] customersSequences;

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

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    // Quantity carried by each truck
    HxExpression[][] quantity;

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

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

        // Sequence of customers visited by each truck
        for (int k = 0; k < nbTrucks; ++k)
            customersSequences[k] = model.List(nbCustomers);

        // Quantity carried by each truck for each customer
        for (int k = 0; k < nbTrucks; ++k)
        {
            quantity[k] = new HxExpression[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
                quantity[k][i] = model.Float(0, demandsData[i]);
        }

        // Every customer must be visited by at least one truck
        model.Constraint(model.Cover(customersSequences));

        // Create HexalyOptimizer arrays to be able to access them with "at" operators
        HxExpression distMatrix = model.Array(distMatrixData);
        HxExpression distDepot = model.Array(distDepotData);

        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 carried in each route must not exceed the truck capacity
            HxExpression quantityArray = model.Array(quantity[k]);
            HxExpression quantityLambda = model.LambdaFunction(j => quantityArray[j]);
            HxExpression routeQuantity = model.Sum(sequence, quantityLambda);
            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(trucksUsed[k], distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);
        }

        for (int i = 0; i < nbCustomers; ++i)
        {
            // Each customer must receive at least its demand
            HxExpression quantityServed = model.Sum();
            for (int k = 0; k < nbTrucks; ++k)
                quantityServed.AddOperand(
                    quantity[k][i] * model.Contains(customersSequences[k], i)
                );
            model.Constraint(quantityServed >= demandsData[i]);
        }

        totalDistance = model.Sum(distRoutes);

        // Objective: minimize the distance traveled
        model.Minimize(totalDistance);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - Each line k contains the customers visited by the truck k */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalDistance.GetValue());
            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
                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: Sdvrp 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 (Sdvrp model = new Sdvrp())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

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

    private void ReadInputSdvrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            splitted = SplitInput(input);

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

            demandsData = new long[nbCustomers];

            splitted = SplitInput(input);
            for (int n = 0; n < nbCustomers; ++n)
                demandsData[n] = int.Parse(splitted[n]);

            splitted = SplitInput(input);
            int depotX = int.Parse(splitted[0]);
            int depotY = int.Parse(splitted[1]);

            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            for (int n = 0; n < nbCustomers; ++n)
            {
                splitted = SplitInput(input);
                customersX[n] = int.Parse(splitted[0]);
                customersY[n] = int.Parse(splitted[1]);
            }

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

    // Compute the distance matrix
    private void ComputeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY)
    {
        distMatrixData = new long[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new long[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                long dist = ComputeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }
        distDepotData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            distDepotData[i] = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
    }

    private static long ComputeDist(int xi, int xj, int yi, int yj)
    {
        double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        return (long)Math.Floor(exactDist + 0.5);
    }
}

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

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    private int truckCapacity;

    // Demand of each customer
    private long[] demandsData;

    // Distance matrix
    private long[][] distMatrixData;

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

    // Number of trucks
    private int nbTrucks;

    // Decision variables
    private HxExpression[] customersSequences;

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

    // Quantity carried by each truck
    private HxExpression[][] quantity;

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

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

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

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

        trucksUsed = new HxExpression[nbTrucks];
        customersSequences = new HxExpression[nbTrucks];
        distRoutes = new HxExpression[nbTrucks];
        quantity = new HxExpression[nbTrucks][nbCustomers];

        // Sequence of customers visited by each truck
        for (int k = 0; k < nbTrucks; ++k)
            customersSequences[k] = model.listVar(nbCustomers);

        // Quantity carried by each truck for each customer
        for (int k = 0; k < nbTrucks; ++k)
            for (int i = 0; i < nbCustomers; ++i)
                quantity[k][i] = model.floatVar(0, demandsData[i]);

        // Every customer must be visited by at least one truck
        model.constraint(model.cover(customersSequences));

        // Create HexalyOptimizer arrays to be able to access them with "at" operators
        HxExpression distMatrix = model.array(distMatrixData);
        HxExpression distDepot = model.array(distDepotData);

        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] = model.gt(c, 0);

            // The quantity carried in each route must not exceed the truck capacity
            HxExpression quantityArray = model.array(quantity[k]);
            HxExpression quantityLambda = model.lambdaFunction(j -> model.at(quantityArray, j));
            HxExpression routeQuantity = model.sum(sequence, quantityLambda);
            model.constraint(model.leq(routeQuantity, truckCapacity));

            // Distance traveled by truck k
            HxExpression distLambda = model
                .lambdaFunction(i -> model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)));
            distRoutes[k] = model.sum(model.sum(model.range(1, c), distLambda),
                model.iif(trucksUsed[k], model.sum(model.at(distDepot, model.at(sequence, 0)),
                    model.at(distDepot, model.at(sequence, model.sub(c, 1)))), 0));
        }

        for (int i = 0; i < nbCustomers; ++i) {
            // Each customer must receive at least its demand
            HxExpression quantityServed = model.sum();
            for (int k = 0; k < nbTrucks; ++k)
                quantityServed.addOperand(model.prod(quantity[k][i], model.contains(customersSequences[k], i)));
            model.constraint(model.geq(quantityServed, demandsData[i]));
        }

        totalDistance = model.sum(distRoutes);

        // Objective: minimize the distance traveled
        model.minimize(totalDistance);

        model.close();

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

        optimizer.solve();
    }

    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            nbCustomers = input.nextInt();
            truckCapacity = input.nextInt();
            nbTrucks = nbCustomers;
            input.nextLine();
            demandsData = new long[nbCustomers];
            for (int n = 1; n <= nbCustomers; ++n) {
                demandsData[n - 1] = input.nextInt();
            }

            input.nextLine();

            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            int depotX = input.nextInt();
            int depotY = input.nextInt();
            input.nextLine();
            for (int n = 1; n <= nbCustomers; ++n) {
                customersX[n - 1] = input.nextInt();
                customersY[n - 1] = input.nextInt();
                input.nextLine();
            }

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

    // Compute the distance matrix
    private void computeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY) {
        distMatrixData = new long[nbCustomers][nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                long dist = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

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

    private static long computeDist(int xi, int xj, int yi, int yj) {
        double exactDist = Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
        return (long) Math.floor(exactDist + 0.5);
    }

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

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

    /* Write the solution in a file with the following format:
     * - Each line k contains the customers visited by the truck k */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(totalDistance.getValue());
            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
                HxCollection customersCollection = customersSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); ++i)
                    output.print((customersCollection.get(i) + 1) + " ");
                output.println();
            }
        }
    }
}