Vehicle Routing with Backhauls (VRPB)

Principles learned

  • Add multiple list decision variables

  • Add a partition constraint

  • Use a lambda expression to compute a sum with a variable number of terms

  • Add ternary conditions

  • Access a multi-dimensional array with an “at” operator

  • Add multiple objectives

Problem

../_images/vrpb.svg

In the Vehicle Routing Problem with Backhauls (VRPB), a fleet of vehicles with uniform capacity housed at the same depot must service customers of two categories: the firsts need to be delivered a known quantity of goods from the single depot while the laters need to send a quantity of goods back to the depot. The model should compute routes representing the order in which customers are visited by a vehicle. Similarly to CVRP, a feasible route must not violate the capacity of a vehicle, and each customer must be served by one route. The specificity of this problem is that all deliveries must be made before any pickups. In a real-life distribution situation where the same vehicles can be used for delivery and pickup, it might be important to introduce this new constraint to ensure the feasability of the goods management. The objective is to minimize the total distance traveled by all vehicles.

Download the example


Data

The instances were provided by LKH-3, and JSprit in the TSPLib format as follows:

  • The number of nodes follows the keyword DIMENSION (there is one warehouse so the number of customers is the number of nodes minus 1).

  • The number of trucks available follows the keyword VEHICLES

  • The truck capacity follows the keyword CAPACITY.

  • The edge type follows EDGE_WEIGHT_TYPE. Note that in our model the only edge type accepted is EXACT_2D.

  • After the keyword NODE_COORD_SECTION, for each node is listed its id and the corresponding x and y coordinates.

  • After the keyword DEMAND_SECTION, for each node is listed its id and the corresponding demand.

  • After the keyword BACKHAUL_SECTION, ids of pickup nodes are listed.

  • Warehouses are listed after the keyword DEPOT_SECTION. Note that in our model only one warehouse is accepted.

Program

The Hexaly model is an extension of the CVRP model. It defines a route for each vehicle as a list variable. Using the partition operator ensures that each customer is assigned to exactly one route.

The precedency constraint is defined in parallel for each vehicles with a for loop. For a specific route, it must ensure that all delivery nodes are visited before any pickup node. A simple way of writing it is to go through each pair of consecutive nodes (i, i+1) of a route, and check that it isn’t a forbidden transition, i.e. node i and node i+1 can’t be respectively a pickup node and a backup node at the same time. This can be written as a basic logic expression using the dictionary backhauls.

Similarly, the distance traveled inside of a route is computed with the operator sum applied with a lambda function to go through all pair of nodes, and access the distance between two nodes precomputed in distanceMatrix.

The capacity constraints are specified using the sum operator on each node of the route by using precomputed arrays deliveryDemand and pickupDemand. Combined with the precedency constraint, they ensure that a feasible solution never exceeds the capacity of a vehicle.

Execution:
hexaly vehicle_routing_backhauls.hxm inFileName=instances/A1.vrpb [hxTimeLimit=] [solFileName=]
use io;

function input() {
    usage = "Usage: hexaly vehicle_routing_backhauls.hxm inFileName=inputFile "
        + "[solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;
    
    readInputVrpb();

    computeDistanceMatrix();
}

function model() {
    // Sequence of customers visited by each truck
    customersSequences[k in 0...nbTrucks] <- list(nbCustomers);
    
    // All customers must be visited by exactly one truck
    constraint partition(customersSequences);

    for [k in 0...nbTrucks] {
        local sequence <- customersSequences[k];
        local c <- count(sequence);

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

        // A pickup cannot be followed by a delivery
        constraint and(1...c, i => isBackhaul[sequence[i-1]] <= isBackhaul[sequence[i]]);

        // The quantity needed in each route must not exceed the truck capacity
        routeDeliveryQuantity <- sum(sequence, j => deliveryDemands[j]);
        constraint routeDeliveryQuantity <= truckCapacity;
        routePickupQuantity <- sum(sequence, j => pickupDemands[j]);
        constraint routePickupQuantity <= truckCapacity;

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

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

    // Total distance traveled
    totalDistance <- sum[k in 0...nbTrucks](routeDistances[k]);

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

/* Parametrize the solver */
function param() { 
    if (hxTimeLimit == nil) hxTimeLimit = 20;
}

/* Write the solution in a file with the following format:
 * - number of trucks used and total distance
 * - for each truck the customers visited (omitting the start/end at the depot) */
function output() {
    if (solFileName == nil) return;
    local outfile = io.openWrite(solFileName);

    outfile.println(nbTrucksUsed.value, " ", totalDistance.value);
    for [k in 0...nbTrucks] {
        if (trucksUsed[k].value != 1) continue;
        // Values in sequence are in 0...nbCustomers.
        // +2 is to put it back in 2...nbCustomers+2
        // as in the data files (1 being the depot)
        for [customer in customersSequences[k].value]
            outfile.print(customer + 2, " ");
        outfile.println();
    }
}

function readInputVrpb() {
    local inFile = io.openRead(inFileName);
    local nbNodes = 0;
    while (true) {
        local str = inFile.readString();
        if (str.startsWith("DIMENSION")) {
            if (!str.endsWith(":")) str = inFile.readString();
            nbNodes = inFile.readInt();
            nbCustomers = nbNodes - 1;
        } else if ((str.startsWith("VEHICLES"))) {
            if (!str.endsWith(":")) str = inFile.readString();
            nbTrucks = inFile.readInt();
        } else if ((str.startsWith("CAPACITY"))) {
            if (!str.endsWith(":")) str = inFile.readString();
            truckCapacity = inFile.readInt();
        } else if ((str.startsWith("EDGE_WEIGHT_TYPE"))) {
            if (!str.endsWith(":")) str = inFile.readString();
            local weightType = inFile.readString();
            if (weightType != "EXACT_2D")
                throw "Edge Weight Type " + weightType + " is not supported (only EXACT_2D)";
        } else if (str.startsWith("NODE_COORD_SECTION")) {
            break;
        } else {
            local dump = inFile.readln();
        }
    }

    // -2 because original customer indices are in 2..nbNodes
    for [n in 0...nbNodes] {
        local node_id = inFile.readInt();
        if (n + 1 != node_id) throw "Unexpected index";
        if (node_id == 1) {
            depotX = round(inFile.readDouble());
            depotY = round(inFile.readDouble());
        } else {
            customersX[node_id - 2] = round(inFile.readDouble());
            customersY[node_id - 2] = round(inFile.readDouble());
        }
    }

    dump = inFile.readln();
    if (!dump.startsWith("DEMAND_SECTION")) throw "Expected keyword DEMAND_SECTION";
    for [n in 1..nbNodes] {
        if (n != inFile.readInt()) throw "Unexpected index";
        local demand = inFile.readInt();
        if (n == 1) {
            if (demand != 0) throw "Demand for depot should be O";
        } else {
            demands[n - 2] = demand; // -2 because original customer indices are in 2..nbNodes
        }
    }
    
    dump = inFile.readln();
    if (!dump.startsWith("BACKHAUL_SECTION")) throw "Expected keyword BACKHAUL_SECTION";
    isBackhaul[i in 0...nbCustomers] = 0; 
    while (true) {
        local node_id = inFile.readInt();
        if (node_id == -1) break;
        // -2 because original customer indices are in 2..nbNodes
        isBackhaul[node_id - 2] = 1;
    }
    for [i in 0...nbCustomers] {
        deliveryDemands[i] = isBackhaul[i] ? 0 : demands[i];
        pickupDemands[i] = isBackhaul[i] ? demands[i] : 0;
    }

    dump = inFile.readln();
    if (!dump.startsWith("DEPOT_SECTION")) throw "Expected keyword DEPOT_SECTION";
    local depotId = inFile.readInt();
    if (depotId != 1) throw "Depot id is supposed to be 1";
    local endOfDepotSection = inFile.readInt();
    if (endOfDepotSection != -1) throw "Expecting only one depot, more than one found";
}

/* Compute the distance matrix */
function computeDistanceMatrix() {
    for [i in 0...nbCustomers] {
        distanceMatrix[i][i] = 0;
        for [j in i+1...nbCustomers] {
            local localDistance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
            distanceMatrix[j][i] = localDistance;
            distanceMatrix[i][j] = localDistance;
        }
    }

    for [i in 0...nbCustomers] {
        local localDistance = computeDist(depotX, customersX[i], depotY, customersY[i]);
        distanceDepot[i] = localDistance;
    }
}

function computeDist(xi, xj, yi, yj) {
    local exactDist = sqrt(pow((xi - xj), 2) + pow((yi - yj), 2));
    return round(exactDist);
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python vehicle_routing_backhauls.py instances\A1.vrpb
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python vehicle_routing_backhauls.py instances/A1.vrpb
import hexaly.optimizer
import sys
import math


def read_elem(filename):
    with open(filename) as f:
        return [str(elem) for elem in f.read().split()]


def main(instance_file, str_time_limit, output_file):

    #
    # Read instance data
    #
    nb_customers, nb_trucks, truck_capacity, dist_matrix_data, dist_depot_data, \
        delivery_demands_data, pickup_demands_data, backhaul_data = read_input_vrpb(
            instance_file)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        model = optimizer.model

        # Sequence of customers visited by each truck
        customers_sequences = [model.list(nb_customers)
                               for _ in range(nb_trucks)]

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

        # Create Hexaly arrays to be able to access them with an "at" operator
        delivery_demands = model.array(delivery_demands_data)
        pickup_demands = model.array(pickup_demands_data)
        dist_matrix = model.array(dist_matrix_data)
        dist_depot = model.array(dist_depot_data)

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

        dist_routes = [None] * nb_trucks
        is_backhaul = model.array(backhaul_data.values())
        for k in range(nb_trucks):
            sequence = customers_sequences[k]
            c = model.count(sequence)

            # A pickup cannot be followed by a delivery
            precedency_lambda = model.lambda_function(lambda i: model.or_(model.not_(
                model.at(is_backhaul, sequence[i-1])), model.at(is_backhaul, sequence[i])))
            model.constraint(model.and_(model.range(1, c), precedency_lambda))

            # The quantity needed in each route must not exceed the truck capacity
            delivery_demand_lambda = model.lambda_function(
                lambda j: delivery_demands[j])
            route_pickup_quantity = model.sum(sequence, delivery_demand_lambda)
            model.constraint(route_pickup_quantity <= truck_capacity)

            pickup_demand_lambda = model.lambda_function(
                lambda j: pickup_demands[j])
            route_pickup_quantity = model.sum(sequence, pickup_demand_lambda)
            model.constraint(route_pickup_quantity <= truck_capacity)

            # Distance traveled by each truck
            dist_lambda = model.lambda_function(lambda i:
                                                model.at(dist_matrix,
                                                         sequence[i - 1],
                                                         sequence[i]))
            dist_routes[k] = model.sum(model.range(1, c), dist_lambda) \
                + model.iif(c > 0,
                            dist_depot[sequence[0]] +
                            dist_depot[sequence[c - 1]],
                            0)

        # Total number of trucks used
        nb_trucks_used = model.sum(trucks_used)

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

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

        model.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = int(str_time_limit)

        optimizer.solve()

        #
        # Write the solution in a file with the following format:
        #  - number of trucks used and total distance
        #  - for each truck the customers visited (omitting the start/end at the depot)
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("%d %d\n" %
                        (nb_trucks_used.value, total_distance.value))
                for k in range(nb_trucks):
                    if trucks_used[k].value != 1:
                        continue
                    # Values in sequence are in 0...nbCustomers. +2 is to put it back
                    # in 2...nbCustomers+2 as in the data files (1 being the depot)
                    for customer in customers_sequences[k].value:
                        f.write("%d " % (customer + 2))
                    f.write("\n")


# The input files follow the "CVRPLib" format
def read_input_vrpb(filename):
    file_it = iter(read_elem(filename))

    nb_nodes = 0
    while True:
        token = next(file_it)
        if token == "DIMENSION":
            next(file_it)  # Removes the ":"
            nb_nodes = int(next(file_it))
            nb_customers = nb_nodes - 1
        elif token == "VEHICLES":
            next(file_it)  # Removes the ":"
            nb_trucks = int(next(file_it))
        elif token == "CAPACITY":
            next(file_it)  # Removes the ":"
            truck_capacity = int(next(file_it))
        elif token == "EDGE_WEIGHT_TYPE":
            next(file_it)  # Removes the ":"
            token = next(file_it)
            if token != "EXACT_2D":
                print("Edge Weight Type " + token +
                      " is not supported (only EXACT_2D)")
                sys.exit(1)
        elif token == "NODE_COORD_SECTION":
            break

    customers_x = [None] * nb_customers
    customers_y = [None] * nb_customers
    depot_x = 0
    depot_y = 0
    for n in range(nb_nodes):
        node_id = int(next(file_it))
        if node_id != n + 1:
            print("Unexpected index")
            sys.exit(1)
        if node_id == 1:
            depot_x = int(next(file_it))
            depot_y = int(next(file_it))
        else:
            # -2 because original customer indices are in 2..nbNodes
            customers_x[node_id - 2] = int(next(file_it))
            customers_y[node_id - 2] = int(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)

    token = next(file_it)
    if token != "DEMAND_SECTION":
        print("Expected token DEMAND_SECTION")
        sys.exit(1)

    demands = [None] * nb_customers
    for n in range(nb_nodes):
        node_id = int(next(file_it))
        if node_id != n + 1:
            print("Unexpected index")
            sys.exit(1)
        if node_id == 1:
            if int(next(file_it)) != 0:
                print("Demand for depot should be 0")
                sys.exit(1)
        else:
            # -2 because original customer indices are in 2..nbNodes
            demands[node_id - 2] = int(next(file_it))

    token = next(file_it)
    if token != "BACKHAUL_SECTION":
        print("Expected token BACKHAUL_SECTION")
        sys.exit(1)

    is_backhaul = {i: False for i in range(nb_customers)}
    while True:
        node_id = int(next(file_it))
        if node_id == -1:
            break
        # -2 because original customer indices are in 2..nbNodes
        is_backhaul[node_id - 2] = True
    delivery_demands = [0 if is_backhaul[i] else demands[i]
                        for i in range(nb_customers)]
    pickup_demands = [demands[i] if is_backhaul[i]
                      else 0 for i in range(nb_customers)]

    token = next(file_it)
    if token != "DEPOT_SECTION":
        print("Expected token DEPOT_SECTION")
        sys.exit(1)

    depot_id = int(next(file_it))
    if depot_id != 1:
        print("Depot id is supposed to be 1")
        sys.exit(1)

    end_of_depot_section = int(next(file_it))
    if end_of_depot_section != -1:
        print("Expecting only one depot, more than one found")
        sys.exit(1)

    return nb_customers, nb_trucks, truck_capacity, distance_matrix, distance_depots, \
        delivery_demands, pickup_demands, is_backhaul


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


# 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 vehicle_routing_backhauls.py input_file [output_file] [time_limit]")
        sys.exit(1)

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

    main(instance_file, str_time_limit, output_file)
Compilation / Execution (Windows)
cl /EHsc vehicle_routing_backhauls.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
vehicle_routing_backhauls instances\A1.vrpb
Compilation / Execution (Linux)
g++ vehicle_routing_backhauls.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o vehicle_routing_backhauls
./vehicle_routing_backhauls instances/A1.vrpb
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each customer
    vector<int> pickupDemandsData;
    vector<int> deliveryDemandsData;

    // Type of each customer
    vector<int> isBackhaulData;

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

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

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    /* Read instance data */
    void readInstance(const string &fileName)
    {
        readInputVrpb(fileName);
    }

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

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

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

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

        trucksUsed.resize(nbTrucks);
        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;

            // A pickup cannot be followed by a delivery
            HxExpression precedencyLambda =
                model.createLambdaFunction([&](HxExpression i)
                                           { return model.leq(model.at(isBackhaul, sequence[i - 1]), model.at(isBackhaul, sequence[i])); });
            model.constraint(model.and_(model.range(1, c), precedencyLambda));

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression deliveryDemandLambda =
                model.createLambdaFunction([&](HxExpression j)
                                           { return deliveryDemands[j]; });
            HxExpression routeDeliveryQuantity = model.sum(sequence, deliveryDemandLambda);
            model.constraint(routeDeliveryQuantity <= truckCapacity);
            HxExpression pickupDemandLambda =
                model.createLambdaFunction([&](HxExpression j)
                                           { return pickupDemands[j]; });
            HxExpression routePickupQuantity = model.sum(sequence, pickupDemandLambda);
            model.constraint(routePickupQuantity <= truckCapacity);

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

        // Total number of trucks used
        nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());

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

        // Objective: minimize the number of trucks used, then minimize the distance traveled
        model.minimize(nbTrucksUsed);
        model.minimize(totalDistance);

        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - number of trucks used and total distance
     *  - for each truck the customers visited (omitting the start/end at the depot) */
    void writeSolution(const string &fileName)
    {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << nbTrucksUsed.getValue() << " " << totalDistance.getValue() << endl;
        for (int k = 0; k < nbTrucks; ++k)
        {
            if (trucksUsed[k].getValue() != 1)
                continue;
            // Values in sequence are in 0...nbCustomers. +2 is to put it back in 2...nbCustomers+2
            // as in the data files (1 being the depot)
            HxCollection customersCollection = customersSequences[k].getCollectionValue();
            for (int i = 0; i < customersCollection.count(); ++i)
            {
                outfile << customersCollection[i] + 2 << " ";
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the "CVRPLib" format
    void readInputVrpb(const string &fileName)
    {
        ifstream infile(fileName.c_str());
        if (!infile.is_open())
        {
            throw std::runtime_error("File cannot be opened.");
        }

        string str;
        char *pch;
        char *line;
        int nbNodes;
        while (true)
        {
            getline(infile, str);
            line = strdup(str.c_str());
            pch = strtok(line, " :");
            if (strcmp(pch, "DIMENSION") == 0)
            {
                pch = strtok(NULL, " :");
                nbNodes = atoi(pch);
                nbCustomers = nbNodes - 1;
            }
            else if (strcmp(pch, "VEHICLES") == 0)
            {
                pch = strtok(NULL, " :");
                nbTrucks = atoi(pch);
            }
            else if (strcmp(pch, "CAPACITY") == 0)
            {
                pch = strtok(NULL, " :");
                truckCapacity = atoi(pch);
            }
            else if (strcmp(pch, "EDGE_WEIGHT_TYPE") == 0)
            {
                pch = strtok(NULL, " :");
                if (strcmp(pch, "EXACT_2D") != 0)
                {
                    throw std::runtime_error("Only Edge Weight Type EXACT_2D is supported");
                }
            }
            else if (strcmp(pch, "NODE_COORD_SECTION") == 0)
            {
                break;
            }
        }

        vector<int> customersX(nbCustomers);
        vector<int> customersY(nbCustomers);
        int depotX, depotY;
        for (int n = 1; n <= nbNodes; ++n)
        {
            int id;
            infile >> id;
            if (id != n)
            {
                throw std::runtime_error("Unexpected index");
            }
            if (n == 1)
            {
                infile >> depotX;
                infile >> depotY;
            }
            else
            {
                // -2 because original customer indices are in 2..nbNodes
                infile >> customersX[n - 2];
                infile >> customersY[n - 2];
            }
        }

        computeDistanceMatrix(depotX, depotY, customersX, customersY);

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " :");
        if (strcmp(pch, "DEMAND_SECTION") != 0)
        {
            throw std::runtime_error("Expected keyword DEMAND_SECTION");
        }

        vector<int> demandsData(nbCustomers);
        for (int n = 1; n <= nbNodes; ++n)
        {
            int id;
            infile >> id;
            if (id != n)
            {
                throw std::runtime_error("Unexpected index");
            }
            int demand;
            infile >> demand;
            if (n == 1)
            {
                if (demand != 0)
                {
                    throw std::runtime_error("Demand for depot should be O");
                }
            }
            else
            {
                // -2 because original customer indices are in 2..nbNodes
                demandsData[n - 2] = demand;
            }
        }

        isBackhaulData.resize(nbCustomers);
        fill(isBackhaulData.begin(), isBackhaulData.end(), 0);
        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " :");
        if (strcmp(pch, "BACKHAUL_SECTION") != 0)
        {
            throw std::runtime_error("Expected keyword BACKHAUL_SECTION");
        }
        while (true)
        {
            int id;
            infile >> id;
            if (id == -1)
                break;
            // -2 because original customer indices are in 2..nbNodes
            isBackhaulData[id - 2] = 1;
        }

        deliveryDemandsData.resize(nbCustomers);
        pickupDemandsData.resize(nbCustomers);
        for (int i = 0; i <= nbCustomers; ++i)
        {
            if (isBackhaulData[i])
            {
                deliveryDemandsData[i] = 0;
                pickupDemandsData[i] = demandsData[i];
            }
            else
            {
                deliveryDemandsData[i] = demandsData[i];
                pickupDemandsData[i] = 0;
            }
        }

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " :");
        if (strcmp(pch, "DEPOT_SECTION") != 0)
        {
            throw std::runtime_error("Expected keyword DEPOT_SECTION");
        }

        int depotId;
        infile >> depotId;
        if (depotId != 1)
        {
            throw std::runtime_error("Depot id is supposed to be 1");
        }

        int endOfDepotSection;
        infile >> endOfDepotSection;
        if (endOfDepotSection != -1)
        {
            throw std::runtime_error("Expecting only one depot, more than one found");
        }

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

    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: vehicle_routing_backhauls 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
    {
        VehicleRoutingWithBackhauls model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    }
    catch (const exception &e)
    {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc VehicleRoutingWithBackhauls.cs /reference:Hexaly.NET.dll
VehicleRoutingWithBackhauls instances\A1.vrpb
using System;
using System.IO;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each customer
    long[] deliveryDemandsData;
    long[] pickupDemandsData;

    // Type of each customer
    int[] isBackhaulData;

    // Distance matrix between customers
    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;

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

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

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputVrpb(fileName);
    }

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

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

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

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

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

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression deliveryDemands = model.Array(deliveryDemandsData);
        HxExpression pickupDemands = model.Array(pickupDemandsData);
        HxExpression isBackhaul = model.Array(isBackhaulData);
        HxExpression distDepot = model.Array(distDepotData);
        HxExpression distMatrix = model.Array(distMatrixData);

        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;

            // A pickup cannot be followed by a delivery
            HxExpression precedencyLambda = model.LambdaFunction(
                i => isBackhaul[sequence[i - 1]] <= isBackhaul[sequence[i]]
            );
            model.Constraint(model.And(model.Range(1, c), precedencyLambda));

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression deliveryDemandLambda = model.LambdaFunction(j => deliveryDemands[j]);
            HxExpression routeDeliveryQuantity = model.Sum(sequence, deliveryDemandLambda);
            model.Constraint(routeDeliveryQuantity <= truckCapacity);
            HxExpression pickupDemandLambda = model.LambdaFunction(j => pickupDemands[j]);
            HxExpression routePickupQuantity = model.Sum(sequence, pickupDemandLambda);
            model.Constraint(routePickupQuantity <= truckCapacity);

            // Distance traveled by truck k
            HxExpression distLambda = model.LambdaFunction(
                i => distMatrix[sequence[i - 1], sequence[i]]
            );
            distRoutes[k] =
                model.Sum(model.Range(1, c), distLambda)
                + model.If(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);
        }

        nbTrucksUsed = model.Sum(trucksUsed);
        totalDistance = model.Sum(distRoutes);

        // Objective: minimize the number of trucks used, then minimize the distance traveled
        model.Minimize(nbTrucksUsed);
        model.Minimize(totalDistance);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - number of trucks used and total distance
     *  - for each truck the customers visited (omitting the start/end at the depot) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(nbTrucksUsed.GetValue() + " " + totalDistance.GetValue());
            for (int k = 0; k < nbTrucks; ++k)
            {
                if (trucksUsed[k].GetValue() != 1)
                    continue;
                // Values in sequence are in 0...nbCustomers. +2 is to put it back in 2...nbCustomers+2
                // as in the data files (1 being the depot)
                HxCollection customersCollection = customersSequences[k].GetCollectionValue();
                for (int i = 0; i < customersCollection.Count(); ++i)
                    output.Write((customersCollection[i] + 2) + " ");
                output.WriteLine();
            }
        }
    }

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

    // The input files follow the "CVRPLib" format
    private void ReadInputVrpb(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            int nbNodes = 0;
            string[] splitted;
            while (true)
            {
                splitted = input.ReadLine().Split(':');
                if (splitted[0].Contains("DIMENSION"))
                {
                    nbNodes = int.Parse(splitted[1]);
                    nbCustomers = nbNodes - 1;
                }
                else if (splitted[0].Contains("VEHICLES"))
                    nbTrucks = int.Parse(splitted[1]);
                else if (splitted[0].Contains("CAPACITY"))
                    truckCapacity = int.Parse(splitted[1]);
                else if (splitted[0].Contains("EDGE_WEIGHT_TYPE"))
                {
                    if (!splitted[1].Trim().Equals("EXACT_2D"))
                        throw new Exception(
                            "Edge Weight Type " + splitted[1] + " is not supported (only EXACT_2D)"
                        );
                }
                else if (splitted[0].Contains("NODE_COORD_SECTION"))
                    break;
            }
            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            int depotX = 0,
                depotY = 0;
            for (int n = 1; n <= nbNodes; ++n)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n)
                    throw new Exception("Unexpected index");
                if (n == 1)
                {
                    depotX = int.Parse(splitted[1]);
                    depotY = int.Parse(splitted[2]);
                }
                else
                {
                    // -2 because original customer indices are in 2..nbNodes
                    customersX[n - 2] = int.Parse(splitted[1]);
                    customersY[n - 2] = int.Parse(splitted[2]);
                }
            }

            ComputeDistanceMatrix(depotX, depotY, customersX, customersY);

            splitted = input.ReadLine().Split(':');
            if (!splitted[0].Contains("DEMAND_SECTION"))
                throw new Exception("Expected keyword DEMAND_SECTION");

            long[] demandsData = new long[nbCustomers];
            for (int n = 1; n <= nbNodes; ++n)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n)
                    throw new Exception("Unexpected index");
                var demand = int.Parse(splitted[1]);
                if (n == 1)
                {
                    if (demand != 0)
                        throw new Exception("Depot demand is supposed to be 0");
                }
                else
                {
                    // -2 because original customer indices are in 2..nbNodes
                    demandsData[n - 2] = demand;
                }
            }

            splitted = input.ReadLine().Split(':');
            if (!splitted[0].Contains("BACKHAUL_SECTION"))
                throw new Exception("Expected keyword BACKHAUL_SECTION");

            isBackhaulData = new int[nbCustomers];
            splitted = input.ReadLine().Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
            foreach (var item in splitted)
            {
                var node_id = int.Parse(item);
                if (node_id == -1)
                    continue;
                // -2 because original customer indices are in 2..nbNodes
                isBackhaulData[node_id - 2] = 1;
            }

            deliveryDemandsData = new long[nbCustomers];
            pickupDemandsData = new long[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
            {
                if (isBackhaulData[i] == 1)
                {
                    pickupDemandsData[i] = demandsData[i];
                }
                else
                {
                    deliveryDemandsData[i] = demandsData[i];
                }
            }

            splitted = input.ReadLine().Split(':');
            if (!splitted[0].Contains("DEPOT_SECTION"))
                throw new Exception("Expected keyword DEPOT_SECTION");

            int depotId = int.Parse(input.ReadLine());
            if (depotId != 1)
                throw new Exception("Depot id is supposed to be 1");

            int endOfDepotSection = int.Parse(input.ReadLine());
            if (endOfDepotSection != -1)
                throw new Exception("Expecting only one depot, more than one found");
        }
    }

    // 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 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 Convert.ToInt64(Math.Round(exactDist));
    }
}
Compilation / Execution (Windows)
javac VehicleRoutingWithBackhauls.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. VehicleRoutingWithBackhauls instances\A1.vrpb
Compilation / Execution (Linux)
javac VehicleRoutingWithBackhauls.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. VehicleRoutingWithBackhauls instances/A1.vrpb
import java.util.*;

import org.w3c.dom.ls.LSException;

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

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

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    private int truckCapacity;

    // Demand on each customer
    private long[] deliveryDemandsData;
    private long[] pickupDemandsData;

    // Type of each customer
    private int[] isBackhaulData;

    // Distance matrix between customers
    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;

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

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

    private VehicleRoutingWithBackhauls(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];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];

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

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

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression deliveryDemands = model.array(deliveryDemandsData);
        HxExpression pickupDemands = model.array(pickupDemandsData);
        HxExpression isBackhaul = model.array(isBackhaulData);
        HxExpression distDepot = model.array(distDepotData);
        HxExpression distMatrix = model.array(distMatrixData);

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

            // A pickup cannot be followed by a delivery
            HxExpression precedencyLambda = model.lambdaFunction(
                    i -> model.leq(
                            model.at(isBackhaul, model.at(sequence, model.sub(i, 1))),
                            model.at(isBackhaul, model.at(sequence, i))));
            model.constraint(model.and(model.range(1, c), precedencyLambda));

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression deliveryDemandLambda = model.lambdaFunction(j -> model.at(deliveryDemands, j));
            HxExpression routeDeliveryQuantity = model.sum(sequence, deliveryDemandLambda);
            model.constraint(model.leq(routeDeliveryQuantity, truckCapacity));
            HxExpression pickupDemandLambda = model.lambdaFunction(j -> model.at(pickupDemands, j));
            HxExpression routePickupQuantity = model.sum(sequence, pickupDemandLambda);
            model.constraint(model.leq(routePickupQuantity, 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(model.gt(c, 0), model.sum(model.at(distDepot, model.at(sequence, 0)),
                            model.at(distDepot, model.at(sequence, model.sub(c, 1)))), 0));
        }

        nbTrucksUsed = model.sum(trucksUsed);
        totalDistance = model.sum(distRoutes);

        // Objective: minimize the number of trucks used, then minimize the distance
        // traveled
        model.minimize(nbTrucksUsed);
        model.minimize(totalDistance);

        model.close();

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

        optimizer.solve();
    }

    /*
     * Write the solution in a file with the following format:
     * - number of trucks used and total distance
     * - for each truck the customers visited (omitting the start/end at the depot)
     */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(nbTrucksUsed.getValue() + " " + totalDistance.getValue());
            for (int k = 0; k < nbTrucks; ++k) {
                if (trucksUsed[k].getValue() != 1)
                    continue;
                // Values in sequence are in 0...nbCustomers. +2 is to put it back in
                // 2...nbCustomers as in the data files (1 being the depot)
                HxCollection customersCollection = customersSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); ++i) {
                    output.print((customersCollection.get(i) + 2) + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "CVRPLib" format
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            int nbNodes = 0;
            String[] splitted;
            while (true) {
                splitted = input.nextLine().split(":");
                if (splitted[0].contains("DIMENSION")) {
                    nbNodes = Integer.parseInt(splitted[1].trim());
                    nbCustomers = nbNodes - 1;
                } else if (splitted[0].contains("VEHICLES")) {
                    nbTrucks = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("CAPACITY")) {
                    truckCapacity = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("EDGE_WEIGHT_TYPE")) {
                    if (splitted[1].trim().compareTo("EXACT_2D") != 0) {
                        throw new RuntimeException(
                                "Edge Weight Type " + splitted[1] + " is not supported (only EXACT_2D)");
                    }
                } else if (splitted[0].contains("NODE_COORD_SECTION")) {
                    break;
                }
            }

            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            int depotX = 0, depotY = 0;
            for (int n = 1; n <= nbNodes; ++n) {
                int id = input.nextInt();
                if (id != n)
                    throw new IOException("Unexpected index");
                if (n == 1) {
                    depotX = input.nextInt();
                    depotY = input.nextInt();
                } else {
                    // -2 because original customer indices are in 2..nbNodes
                    customersX[n - 2] = input.nextInt();
                    customersY[n - 2] = input.nextInt();
                }
            }

            computeDistanceMatrix(depotX, depotY, customersX, customersY);

            input.nextLine().split(":"); // End the last line
            splitted = input.nextLine().split(":");
            if (!splitted[0].contains("DEMAND_SECTION")) {
                throw new RuntimeException("Expected keyword DEMAND_SECTION");
            }

            long[] demandsData = new long[nbCustomers];
            for (int n = 1; n <= nbNodes; ++n) {
                int id = input.nextInt();
                if (id != n)
                    throw new IOException("Unexpected index");
                int demand = input.nextInt();
                if (n == 1) {
                    if (demand != 0)
                        throw new IOException("Depot demand is supposed to be 0");
                } else {
                    // -2 because original customer indices are in 2..nbNodes
                    demandsData[n - 2] = demand;
                }
            }

            input.nextLine().split(":"); // End the last line
            splitted = input.nextLine().split(":");
            if (!splitted[0].contains("BACKHAUL_SECTION")) {
                throw new RuntimeException("Expected keyword DEPOT_SECTION");
            }

            isBackhaulData = new int[nbCustomers];
            while (true) {
                int id = input.nextInt();
                if (id == -1)
                    break;
                // -2 because original customer indices are in 2..nbNodes
                isBackhaulData[id - 2] = 1;
            }

            deliveryDemandsData = new long[nbCustomers];
            pickupDemandsData = new long[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i) {
                if (isBackhaulData[i] == 1) {
                    pickupDemandsData[i] = demandsData[i];
                } else {
                    deliveryDemandsData[i] = demandsData[i];
                }
            }

            input.nextLine().split(":"); // End the last line
            splitted = input.nextLine().split(":");
            if (!splitted[0].contains("DEPOT_SECTION")) {
                throw new RuntimeException("Expected keyword DEPOT_SECTION");
            }

            int depotId = input.nextInt();
            if (depotId != 1)
                throw new IOException("Depot id is supposed to be 1");

            int endOfDepotSection = input.nextInt();
            if (endOfDepotSection != -1) {
                throw new RuntimeException("Expecting only one depot, more than one found");
            }
        }
    }

    // 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 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 Math.round(exactDist);
    }

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

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