Multi-Trip Capacitated Vehicle Routing Problem (MTCVRP)

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):
    nb_customers, nb_trucks, truck_capacity, dist_matrix_data, nb_depots, \
    nb_depot_copies, nb_total_locations, demands_data, max_dist = read_input_multi_trip_vrp(instance_file)
    
    with hexaly.optimizer.HexalyOptimizer() as optimizer:

        #
        # Declare the optimization model
        #
        model = optimizer.model

        # Locations visited by each truck (customer or depot)
        # Add copies of the depots (so that they can be visited multiple times)
        # Add an extra fictive truck (who will visit every depot that will not be visited by real trucks)
        visit_orders = [model.list(nb_total_locations) for k in range(nb_trucks+1)]

        # The fictive truck cannot visit customers
        for i in range(nb_customers):
            model.constraint((model.contains(visit_orders[nb_trucks],i)) == False)

        # All customers must be visited by exactly one truck
        model.constraint(model.partition(visit_orders))

        # Create Hexaly arrays to be able to access them with an "at" operator
        demands = model.array(demands_data)
        dist_matrix = model.array(dist_matrix_data)
        
         # A truck is used if it visits at least one customer
        trucks_used = [(model.count(visit_orders[k]) > 0) for k in range(nb_trucks)]

        dist_routes = [None] * nb_trucks
        for k in range(nb_trucks):
            sequence = visit_orders[k]
            c = model.count(sequence)

            # Compute the quantity in the truck at each step
            route_quantity_lambda = model.lambda_function(lambda i,prev: \
                model.iif(sequence[i] < nb_customers, prev+demands[sequence[i]],0))
            route_quantity = model.array(model.range(0, c), route_quantity_lambda, 0)

            # Trucks cannot carry more than their capacity
            quantity_lambda = model.lambda_function(
                lambda i: route_quantity[i] <= truck_capacity)
            model.constraint(model.and_(model.range(0, c), quantity_lambda))

            # 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,
                            model.at(dist_matrix,nb_customers,sequence[0]) +\
                            model.at(dist_matrix,sequence[c-1],nb_customers),\
                            0)
            
            # Limit distance traveled
            model.constraint( dist_routes[k] <= max_dist)

        # Total distance traveled
        total_distance = model.sum(dist_routes)

        # 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 solution output
        if output_file != None:
            with open(output_file, 'w') as file:
                file.write("File name: %s; totalDistance = %d \n" % (instance_file,total_distance.value))
                for k in range(nb_trucks):
                    if trucks_used[k].value:
                        file.write("Truck %d : " % (k))
                        for customer in visit_orders[k].value:
                            file.write( "%d" % (customer) if customer<nb_customers else "%d" % (-(math.floor((customer-nb_customers)/nb_depot_copies) + 1)))
                            file.write(" ")
                        file.write("\n")


def read_input_multi_trip_vrp(filename):
    if filename.endswith(".dat"):
        return read_input_multi_trip_vrp_dat(filename)
    else:
        raise Exception("Unknown file format")

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

    nb_customers = int(next(file_it))
    nb_depots = int(next(file_it))

    depots_x = [None] * nb_depots
    depots_y = [None] * nb_depots
    for i in range(nb_depots):
        depots_x[i] = int(next(file_it))
        depots_y[i] = int(next(file_it))

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

    truck_capacity = int(next(file_it))//2

    # Skip depots capacity infos (not related to the problem)
    for i in range(nb_depots):
        next(file_it)

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

    nb_depot_copies = 20

    nb_total_locations  = nb_customers + nb_depots*nb_depot_copies

    max_dist = 400

    nb_trucks = 3

    dist_matrix_data = compute_distance_matrix(depots_x, depots_y, customers_x, customers_y, nb_depot_copies)

    return  nb_customers, nb_trucks, truck_capacity, dist_matrix_data, nb_depots, \
        nb_depot_copies, nb_total_locations, demands_data, max_dist


# Compute the distance matrix
def compute_distance_matrix(depots_x, depots_y, customers_x, customers_y, nb_depot_copies):
    nb_customers = len(customers_x)
    nb_depots = len(depots_x)
    nb_total_locations = nb_customers + nb_depots*nb_depot_copies
    dist_matrix = [[0 for _ in range(nb_total_locations)] for _ in range(nb_total_locations)]
    for i in range(nb_customers):
        dist_matrix[i][i] = 0
        for j in range(i,nb_customers):
            dist = compute_dist(customers_x[i], customers_x[j], customers_y[i], customers_y[j])
            dist_matrix[i][j] = dist
            dist_matrix[j][i] = dist
        for d in range(nb_depots):
            dist = compute_dist(customers_x[i], depots_x[d], customers_y[i], depots_y[d])
            for c in range(nb_depot_copies):
                j = nb_customers+d*nb_depot_copies + c
                dist_matrix[i][j] = dist
                dist_matrix[j][i] = dist
    
    for i in range(nb_customers, nb_total_locations):
        for j in range(nb_customers, nb_total_locations):
            # Going from one depot to an other is never worth it
            dist_matrix[i][j] = 100000

    return dist_matrix


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 cvrp.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 MultiTripVRP {
public:
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Number of depots
    int nbDepots;

    // Number of depot copies
    int nbDepotCopies;

    // Total number of locations (customers, depots, depots copies)
    int nbTotalLocations;

    // Capacity of the trucks
    int truckCapacity;

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

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

    // Number of trucks
    int nbTrucks;

    // Maximum distance traveled by a truck
    int maxDist;

    // Decision variables
    vector<HxExpression> visitOrders;

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    void readInstance(const char* fileName) { readInputMultiTripVRP(fileName); }

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

        // Locations visited by each truck (Customers and Depots)
        // Add copies of the depots (so that they can be visited multiple times)
        // Add an extra fictive truck (who will visit every depot that will not be visited by real trucks)
        visitOrders.resize(nbTrucks + 1);
        for (int k = 0; k < nbTrucks + 1; ++k) {
            visitOrders[k] = model.listVar(nbTotalLocations);
        }

        // The fictive truck cannot visit customers
        for(int i=0; i < nbCustomers; i++){
            model.constraint(!(model.contains(visitOrders[nbTrucks],i)));
        }

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

        // Create Hexaly arrays to be able to access them with an "at" operator
        HxExpression demands = model.array(demandsData.begin(), demandsData.end());
        HxExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers + nbDepots * nbDepotCopies; ++n) {
            distMatrix.addOperand(model.array(distMatrixData[n].begin(), distMatrixData[n].end()));
        }

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

        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = visitOrders[k];
            HxExpression c = model.count(sequence);

            // A truck is used if it visits at least one customer
            trucksUsed[k] = c > 0;

            // Compute the quantity in the truck at each step
            HxExpression routeQuantityLambda =
                model.createLambdaFunction([&](HxExpression i, HxExpression prev) { return model.iif(sequence[i] < nbCustomers, prev+demands[sequence[i]],0); });
            HxExpression routeQuantity = model.array(model.range(0, c), routeQuantityLambda, 0);
            // Trucks cannot carry more than their capacity
            HxExpression quantityLambda = model.createLambdaFunction([&](HxExpression i){ return routeQuantity[i]<=truckCapacity; });
            model.constraint(model.and_(model.range(0,c), quantityLambda));

            // 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,
                                model.at(distMatrix,nbCustomers,sequence[0]) +
                                model.at(distMatrix,sequence[c-1],nbCustomers),
                                0);
            
            // Limit distance traveled
            model.constraint(distRoutes[k] <= maxDist);
        }

        // Total distance traveled
        totalDistance = model.sum(distRoutes.begin(), distRoutes.end());

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

        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file */
    void writeSolution(const char* inFile, const string& solFile) {
        ofstream file;
        file.exceptions(ofstream::failbit | ofstream::badbit);
        file.open(solFile.c_str());
        file << "File name: " << inFile << "; total distance = " << totalDistance.getValue() << endl;
        for (int r = 0; r < nbTrucks; ++r) {
            if (trucksUsed[r].getValue()) {
                HxCollection visitCollection = visitOrders[r].getCollectionValue();
                    file << "Truck " << r << " : " ;
                for (hxint i = 0; i < visitCollection.count(); ++i) {
                    file << (visitCollection[i]<nbCustomers ? visitCollection[i] : -(floor((visitCollection[i]-nbCustomers)/nbDepotCopies) + 1)) << " ";
                }
                file << endl;
            }
        }
    }

private:

    void readInputMultiTripVRP(const char* fileName) {
        vector<int> customersX, customersY, depotsX, depotsY;
        string file = fileName;
        ifstream infile(file.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }
        infile >> nbCustomers;
        customersX.resize(nbCustomers);
        customersY.resize(nbCustomers);
        demandsData.resize(nbCustomers);
        distMatrixData.resize(nbTotalLocations);
        infile >> nbDepots;
        depotsX.resize(nbDepots);
        depotsY.resize(nbDepots);
        for (int i = 0; i < nbDepots; ++i) {
            infile >> depotsX[i];
            infile >> depotsY[i];
        }
        for (int i = 0; i < nbCustomers; ++i) {
            infile >> customersX[i];
            infile >> customersY[i];
        }
        infile >> truckCapacity;
        truckCapacity /= 2;
        // Ignore depot infos
        int temp;
        for (int i = 0; i < nbDepots; ++i) {
            infile>>temp;
        }
        for (int i = 0; i < nbCustomers; ++i) {
            infile >> demandsData[i];
        }

        nbDepotCopies = 20;

        nbTotalLocations = nbCustomers + nbDepots*nbDepotCopies;

        maxDist = 400;
            
        nbTrucks = 3;

        computeDistanceMatrix(depotsX, depotsY, customersX, customersY);
    }

    // Compute the distance matrix
    void computeDistanceMatrix(const vector<int>& depotsX, const vector<int>& depotsY, const vector<int>& customersX, const vector<int>& customersY) {
        distMatrixData.resize(nbTotalLocations);
        for (int i = 0; i < nbTotalLocations; ++i) {
            distMatrixData[i].resize(nbTotalLocations);
        }
        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;
            }

            for (int d=0; d<nbDepots; ++d){
                int distance = computeDist(customersX[i], depotsX[d], customersY[i], depotsY[d]);
                for (int c=0; c<nbDepotCopies; ++c){
                    int j = nbCustomers + d*nbDepotCopies + c;
                    distMatrixData[i][j] = distance;
                    distMatrixData[j][i] = distance;
                }
            }
        }

        for(int i=nbCustomers; i<nbTotalLocations; i++){
            for(int j=nbCustomers; j<nbTotalLocations;j++){
                // Going from one depot to an other is never worth it
                distMatrixData[i][j] = 100000;
            }
        }
    }

    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: ./lrp 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 {
        MultiTripVRP model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(instanceFile, solFile);
    } catch (const std::exception& e) {
        std::cerr << "An error occured: " << e.what() << endl;
    }

    return 0;
}

using System;
using System.IO;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

    // Number of depots
    int nbDepots;

    // Number of depot copies
    int nbDepotCopies;

    // Total number of locations (customers, depots, depots copies)
    int nbTotalLocations;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each customer
    int[] demandsData;

    // Distance matrix between customers
    int[][] distMatrixData;

    // Number of trucks
    int nbTrucks;

    // Maximum distance traveled by a truck
    int maxDist;

    // Decision variables
    HxExpression[] visitOrders;

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

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

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

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            nbDepotCopies = 20;
            maxDist = 400;
            string[] line;
            line = ReadNextLine(input);
            nbCustomers = int.Parse(line[0]);
            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            demandsData = new int[nbCustomers];

            line = ReadNextLine(input);
            nbDepots = int.Parse(line[0]);
            int[] depotsX = new int[nbDepots];
            int[] depotsY = new int[nbDepots];

            line = ReadNextLine(input);
            for (int i = 0; i < nbDepots; ++i)
            {
                line = ReadNextLine(input);
                depotsX[i] = int.Parse(line[0]);
                depotsY[i] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);

            for (int i = 0; i < nbCustomers; ++i)
            {
                line = ReadNextLine(input);
                customersX[i] = int.Parse(line[0]);
                customersY[i] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            truckCapacity = int.Parse(line[0]) / 2;
            line = ReadNextLine(input);

            // Ignore depot infos
            for (int i = 0; i < nbDepots; ++i)
            {
                ReadNextLine(input);
            }
            line = ReadNextLine(input);
            for (int i = 0; i < nbCustomers; ++i)
            {
                line = ReadNextLine(input);
                demandsData[i] = int.Parse(line[0]);
            }
            line = ReadNextLine(input);

            nbTotalLocations = nbCustomers + nbDepots * nbDepotCopies;

            nbTrucks = 3;

            ComputeDistanceMatrix(depotsX, depotsY, customersX, customersY);
        }
    }

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

        trucksUsed = new HxExpression[nbTrucks];
        visitOrders = new HxExpression[nbTrucks + 1];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];

        // Locations visited by each truck (Customers and Depots)
        // Add copies of the depots (so that they can be visited multiple times)
        // Add an extra fictive truck (who will visit every depot that will not be visited by real trucks)
        for (int k = 0; k < nbTrucks + 1; ++k)
            visitOrders[k] = model.List(nbTotalLocations);

        // The fictive truck cannot visit customers
        for (int i = 0; i < nbCustomers; i++)
            model.Constraint(!(model.Contains(visitOrders[nbTrucks], i)));

        // All customers must be visited by exactly one truck
        model.Constraint(model.Partition(visitOrders));

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.Array(demandsData);
        HxExpression distMatrix = model.Array(distMatrixData);

        for (int k = 0; k < nbTrucks; ++k)
        {
            HxExpression sequence = visitOrders[k];
            HxExpression c = model.Count(sequence);

            // A truck is used if it visits at least one customer
            trucksUsed[k] = c > 0;

            // Compute the quantity in the truck at each step
            HxExpression routeQuantityLambda = model.LambdaFunction((i, prev) =>
                model.If(sequence[i] < nbCustomers, prev + demands[sequence[i]], 0));
            HxExpression routeQuantity = model.Array(model.Range(0, c), routeQuantityLambda, 0);
            // Trucks cannot carry more than their capacity
            HxExpression quantityLambda = model.LambdaFunction(i => routeQuantity[i] <= truckCapacity);
            model.Constraint(model.And(model.Range(0, c), quantityLambda));

            // 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, distMatrix[nbCustomers, sequence[0]] + distMatrix[sequence[c - 1], nbCustomers], 0);

            model.Constraint(distRoutes[k] <= maxDist);
        }

        // Total distance traveled
        totalDistance = model.Sum(distRoutes);

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

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file */
    void WriteSolution(string infile, string outfile)
    {
        using (StreamWriter output = new StreamWriter(outfile))
        {
            output.WriteLine(
                "File name: " + infile + "; total distance = " + totalDistance.GetValue()
            );
            for (int r = 0; r < nbTrucks; ++r)
            {
                if (trucksUsed[r].GetValue() == 1)
                {
                    output.Write(
                        "Truck " + r + ": "
                    );
                    HxCollection visitCollection = visitOrders[r].GetCollectionValue();
                    for (int i = 0; i < visitCollection.Count(); ++i)
                        output.Write(
                            (visitCollection[i] < nbCustomers ? visitCollection[i] : -(Math.Floor(Convert.ToDouble((visitCollection[i] - nbCustomers) / nbDepotCopies) + 1)))
                            + " ");
                    output.WriteLine();
                }
            }
        }
    }

    String[] ReadNextLine(StreamReader input)
    {
        return input.ReadLine().Split(new[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
    }


    private int computeDist(int xi, int xj, int yi, int yj)
    {
        double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        return Convert.ToInt32(Math.Round(exactDist));
    }

    // Compute the distance matrix
    private void ComputeDistanceMatrix(int[] depotsX, int[] depotsY, int[] customersX, int[] customersY)
    {
        distMatrixData = new int[nbTotalLocations][];
        for (int i = 0; i < nbCustomers + nbDepots * nbDepotCopies; ++i)
            distMatrixData[i] = new int[nbTotalLocations];
        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;
            }

            for (int d = 0; d < nbDepots; d++)
            {
                int distance = computeDist(customersX[i], depotsX[d], customersY[i], depotsY[d]);
                for (int c = 0; c < nbDepotCopies; c++)
                {
                    int j = nbCustomers + d * nbDepotCopies + c;
                    distMatrixData[i][j] = distance;
                    distMatrixData[j][i] = distance;
                }
            }
        }

        for (int i = nbCustomers; i < nbTotalLocations; i++)
        {
            // Going from one depot to an other is never worth it
            for (int j = nbCustomers; j < nbTotalLocations; j++)
                distMatrixData[i][j] = 100000;
        }
    }

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

};

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

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

    // Number of customers
    private int nbCustomers;

    // Number of depots
    private int nbDepots;

    // Number of depot copies
    private int nbDepotCopies;

    // Total number of locations (customers, depots, depots copies)
    private int nbTotalLocations;

    // Capacity of the trucks
    private int truckCapacity;

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

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

    // Number of trucks
    private int nbTrucks;

    // Maximum distance traveled by a truck
    private int maxDist;

    // Decision variables
    private HxExpression[] visitOrders;

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

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

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

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

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

        trucksUsed = new HxExpression[nbTrucks];
        visitOrders = new HxExpression[nbTrucks + 1];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];

        // Locations visited by each truck (Customers and Depots)
        // Add copies of the depots (so that they can be visited multiple times)
        // Add an extra fictive truck (who will visit every depot/ depot copy that will
        // not be visited by real trucks)
        for (int k = 0; k < nbTrucks + 1; ++k)
            visitOrders[k] = m.listVar(nbTotalLocations);

        // The fictive truck cannot visit customers
        for (int k = 0; k < nbCustomers; ++k)
            m.constraint(m.not((m.contains(visitOrders[nbTrucks], k))));

        // All customers must be visited by exactly one truck
        m.constraint(m.partition(visitOrders));

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = m.array(demandsData);
        HxExpression distMatrix = m.array(distMatrixData);

        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = visitOrders[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 routeQuantityLambda = m.lambdaFunction((i, prev) -> m
                    .iif(m.leq(m.at(sequence, i), nbCustomers - 1), m.sum(prev, m.at(demands, m.at(sequence, i))), 0));
            HxExpression routeQuantity = m.array(m.range(0, c), routeQuantityLambda, 0);
            // Trucks cannot carry more than their capacity
            HxExpression quantityLambda = m.lambdaFunction(i -> m.leq(m.at(routeQuantity, i), truckCapacity));
            m.constraint(m.and(m.range(0, c), quantityLambda));

            // 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(m.at(distMatrix, nbCustomers), m.at(sequence, 0)),
                            m.at(m.at(distMatrix, m.at(sequence, m.sub(c, 1))), nbCustomers)), 0));

            m.constraint(m.leq(distRoutes[k], maxDist));
        }

        totalDistance = m.sum(distRoutes);

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

        m.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file */
    private void writeSolution(String infile, String outfile) throws IOException {
        try (PrintWriter output = new PrintWriter(outfile)) {
            output.println("File name: " + infile + "; total distance = " + totalDistance.getValue());
            for (int r = 0; r < nbTrucks; ++r) {
                if (trucksUsed[r].getValue() != 0) {
                    output.print("Truck " + r + " : ");
                    HxCollection visitCollection = visitOrders[r].getCollectionValue();
                    for (int i = 0; i < visitCollection.count(); ++i) {
                        output.print(
                                (visitCollection.get(i) < nbCustomers ? visitCollection.get(i)
                                        : -((int) (Math.floor((visitCollection.get(i) - nbCustomers) / nbDepotCopies)
                                                + 1)))
                                        + " ");
                    }
                    output.println();
                }
            }
        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    private void readInstanceMultiTripVrp(String fileName) throws IOException {
        try (Scanner infile = new Scanner(new File(fileName))) {
            nbDepotCopies = 20;
            nbCustomers = infile.nextInt();
            long[] customersX = new long[nbCustomers];
            long[] customersY = new long[nbCustomers];
            demandsData = new long[nbCustomers];
            nbDepots = infile.nextInt();
            long[] depotsX = new long[nbDepots];
            long[] depotsY = new long[nbDepots];

            for (int i = 0; i < nbDepots; ++i) {
                depotsX[i] = infile.nextInt();
                depotsY[i] = infile.nextInt();
            }
            for (int i = 0; i < nbCustomers; ++i) {
                customersX[i] = infile.nextInt();
                customersY[i] = infile.nextInt();
            }
            truckCapacity = infile.nextInt() / 2;
            for (int i = 0; i < nbDepots; ++i) {
                infile.next();
            }
            for (int i = 0; i < nbCustomers; ++i) {
                demandsData[i] = infile.nextInt();
            }

            nbTotalLocations = nbCustomers + nbDepots * nbDepotCopies;

            nbTrucks = 3;

            maxDist = 400;

            computeDistanceMatrix(depotsX, depotsY, customersX, customersY);

        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    // Compute the distance matrix
    private void computeDistanceMatrix(long[] depotsX, long[] depotsY, long[] customersX, long[] customersY) {
        distMatrixData = new long[nbTotalLocations][nbTotalLocations];

        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                long distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = distance;
                distMatrixData[j][i] = distance;
            }

            for (int d = 0; d < nbDepots; d++) {
                long distance = computeDist(customersX[i], depotsX[d], customersY[i], depotsY[d]);
                for (int c = 0; c < nbDepotCopies; c++) {
                    int j = nbCustomers + d * nbDepotCopies + c;
                    distMatrixData[i][j] = distance;
                    distMatrixData[j][i] = distance;
                }
            }
        }

        for (int i = nbCustomers; i < nbTotalLocations; i++) {
            for (int j = nbCustomers; j < nbTotalLocations; j++) {
                // Going from one depot to an other is never worth it
                distMatrixData[i][j] = 100000;
            }
        }
    }

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

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

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

            MultiTripVrp model = new MultiTripVrp(optimizer);
            model.readInstanceMultiTripVrp(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (strOutfile != null)
                model.writeSolution(instanceFile, strOutfile);
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}