This page is for an old version of Hexaly Optimizer. We recommend that you update your version and read the documentation for the latest stable release.

Vehicule routing problemΒΆ

Principles learnedΒΆ

  • Add multiple list decision variables
  • Add a partition contraint
  • Handle the access of a list after its last element
  • Define a sequence of expressions
  • Access a multi-dimensional array with an “at” operator
  • Add multiple objectives

ProblemΒΆ

../_images/vrp.png

In the capacitated vehicle routing problem (CVRP), a fleet of delivery vehicules with uniform capacity must service customers with known demand for a single commodity. The vehicules start and end their routes at a common depot. Each customer can only be served by one vehicule. The objectives are to minimize the fleet size and assign a sequence of customers to each truck of the fleet minimizing the total distance travelled such that all customers are served and the total demand served by each trcuk does not exceed its capacity.

Download the example

DataΒΆ

The instances provided come from the Augerat et al. Set A instances. They follow the TSPLib format.

The format of the data files is 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 truck capacity follows the keyword CAPACITY.
  • The edge type follows EDGE_WEIGHT_TYPE. Note that in our model the only edge type accepted is EUC_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.
  • Warehouses are listed after the keyword DEPOT_SECTION. Note that in our model only one warehouse is accepted.

The number of available trucks can be defined through the command line. If not, it is deduced from the instance file name that follow this pattern: A-nXX- kNBTRUCKS.vrp

ProgramΒΆ

This LocalSolver model defines a route for each truck k as the list variable customersSequences[k]. It corresponds to the sequence of customers visited. To ensure that all customers must be served, all the list variables are constrained to form a partition, thanks to the “partition” operator.

Along each tour k, the nth node visited is computed by nodeOnVisits[k][n]. By convention, the depot node has the id 0, and customer nodes have ids in [1..nbCustomers]. Note that the list variables customersSequences take their values in [0..nbCustomers-1]. The computation is as follows:

  • The first (n=0) and last (n=nbCustomers+1) visited nodes are the depot: nodeOnVisits[k][n] = 0
  • When n is lower than the number of customers in the sequence, the node visited is the customer id: nodeOnVisits[k][n] <- customersSequences[k][n-1]+1 (We get the index n-1 of customersSequences because n is in [1..nbCustomers] whereas the list can be accessed with the “at” operator on indices in [0..nbCustomers-1]. We then add 1 in order to bring the range of list variables from [0..nbCustomers-1] to [1..nbCustomers] to respect the ids convention).
  • When n is greater or equal than the number of customers in the sequence, the node visited is the depot because the tour is over. Because n exceeds the number of elements in the list, customersSequence[k][n] will return -1 as stipulated in the list variables documentation. Thus, the expression is also nodeOnVisits[k][n] <- customersSequences[k][n-1]+1.

The number of trucks used for the fleet is defined by the number of trucks that serve at least one customer (if their list variable has at least one element). The definition of these expressions is really straightforward thanks to the “count” and the “greater than” operators.

The quantity delivered on each visit is the demand on the node of this visit. This expression is just defined with an “at” operator to access the array demands on the index nodeOnVisits[k][n]. The total quantity delivered by each truck is constrained to be lower than its capacity.

For each truck, the distance travelled from the nth visit to the visit n+1 is defined by distanceNodes[k][n]. It is simply done by applying an “at” operator to the multi-dimensional array distanceMatrix, with the first index nodeOnVisits[k][n] and the second index nodeOnVisits[k][n+1].

Both objectives are defined in lexicographical order: we first minimize the number of trucks used, and then we minimize the total distance travelled by all the trucks.

At the end of this page, below the models, are some insights to adapt the CVRP model to a model solving a CVRPTW.

Execution:
localsolver cvrp.lsp inFileName=instances/A-n32-k5.vrp [lsTimeLimit=] [solFileName=]
/********** cvrp.lsp **********/

use io;

/* Reads instance data
 * The input files follow the "Augerat" format.*/
function input() {
    usage = "\nUsage: localsolver cvrp.lsp " + "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit] [nbTrucks=value]\n";

    if (inFileName == nil) throw usage;
    
    readInputCvrp();
    
    // The number of trucks is usually given in the name of the file
    // nbTrucks can also be given in command line
    if(nbTrucks == nil) nbTrucks = getNbTrucks();

    // Compute distance matrix
    computeDistanceMatrix();
}
 
/* Declares the optimization model. */   
function model() {
    // Sequence of customers visited by each truck.
    // When n >= count(customersSequences[k]) (after the last visited
    // customer), customersSequences[k][n] = -1
    customersSequences[k in 1..nbTrucks] <- list(nbCustomers);

    // All customers must be visited by  the trucks
    constraint partition[k in 1..nbTrucks](customersSequences[k]);

    // Each truck starts (n=0) and ends (n=nbCustomers+1) at the depot (index 0)
    nodeOnVisits[k in 1..nbTrucks][0] = 0;
    nodeOnVisits[k in 1..nbTrucks][nbCustomers+1] = 0;

    // During the route, the actual node visited is 1 + customersSequences[k][n]:
    //  - When n < count(customersSequences[k]) (when a customer is on the nth 
    //      position), nodeOnVisits[k][n+1] = cutomerId. (actual customers ids are 
    //      in [1..nbCustomers], 0 being the warehouse)
    //  - When n >= count(customersSequences[k]) (when the last customer has already 
    //      been visited), nodeOnVisits[k][n+1] = 0 (the warehouse id)
    nodeOnVisits[k in 1..nbTrucks][n in 1..nbCustomers] <- 1 + customersSequences[k][n-1];

    // A truck is used if it visits at least one customer
    trucksUsed[k in 1..nbTrucks] <- count(customersSequences[k]) > 0;
    nbTrucksUsed <- sum[k in 1..nbTrucks](trucksUsed[k]);

    // The quantity needed in each route must not exceed the truck capacity
    for [k in 1..nbTrucks] {
        routeQuantity <- sum[n in 1..nbCustomers](demands[nodeOnVisits[k][n]]);
        constraint routeQuantity <= truckCapacity;
    }

    // Distance traveled from node n to node n+1 by truck k
    distanceNodes[k in 1..nbTrucks][n in 0..nbCustomers] <- distanceMatrix[nodeOnVisits[k][n]][nodeOnVisits[k][n+1]];

    // Distance traveled by each truck
    routeDistances[k in 1..nbTrucks] <- sum[n in 0..nbCustomers](distanceNodes[k][n]);

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

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

/* Parameterizes the solver. */
function param() {
    if(lsTimeLimit == nil) lsTimeLimit=20;
}

/* Writes the solution in a file with the following format :
 *  - number of trucks used and total distance
 *  - for each truck the nodes 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 1..nbTrucks] {
        if (trucksUsed[k].value != 1) continue;
        // Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in [2..nbCustomers+1]
        // as in the data files (1 being the depot)
        for [customer in customersSequences[k].value]
            outfile.print(customer + 2, " ");
        outfile.println();
    }
}


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

    for[n in 1..nbNodes] {
        if (n != inFile.readInt()) throw "Unexpected index";
        nodesX[n-1] = round(inFile.readDouble());
        nodesY[n-1] = 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";
        // demands must start at 0 to be accessed by an "at" operator. Thus
        // node ids will start at 0 in the model.
        demands[n-1] = inFile.readInt();
    }

    dump = inFile.readln();
    if (!dump.startsWith("DEPOT_SECTION")) throw "Expected keyword DEPOT_SECTION";
    local warehouseId = inFile.readInt();
    if (warehouseId != 1) throw "Warehouse id is supposed to be 1";
    local endOfDepotSection = inFile.readInt();
    if (endOfDepotSection != -1) throw "Expecting only one warehouse, more than one found";
    if (demands[0] != 0) throw "Warehouse demand is supposed to be 0";
}

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

function computeDist(i,j) {
    local exactDist = sqrt(pow((nodesX[i] - nodesX[j]), 2) + pow((nodesY[i] - nodesY[j]), 2));
    return round(exactDist);
}

function getNbTrucks() {
    local splitted = split(inFileName,"-k");
    if (count(splitted) >= 2) {
        local numvrp = splitted[count(splitted)-1];
        splitted = split(numvrp, ".");
        if (count(splitted) == 2) return toInt(splitted[0]);
    } else {
        println("Error: nbTrucks could not be read from the file name. Enter it from the command line");
        throw usage;
    }
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python27\
python cvrp.py instances\A-n32-k5.vrp
Execution (Linux)
export PYTHONPATH=/opt/localsolver_XXX/bin/python27/
python cvrp.py instances/A-n32-k5.vrp
########## cvrp.py ##########

import localsolver
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, str_nb_trucks):
    nb_trucks = int(str_nb_trucks)

    #
    # Reads instance data
    #
    (nb_customers, truck_capacity, distance_matrix, demands) = read_input_cvrp(instance_file)
    
    # The number of trucks is usually given in the name of the file
    # nb_trucks can also be given in command line
    if nb_trucks == 0:
        nb_trucks = get_nb_trucks(instance_file)

    with localsolver.LocalSolver() as ls:
        #
        # Declares the optimization model
        #
        model = ls.model

        # Sequence of customers visited by each truck.
        # When n >= count(customers_sequences[k]) (after the last visited
        # customer), customers_sequences[k][n] = -1
        customers_sequences = [model.list(nb_customers) for k in range(nb_trucks)]

        # All customers must be visited by  the trucks
        model.constraint(model.partition(customers_sequences))

        # Each truck starts (n=0) and ends (n=nb_customers+1) at the depot (index 0)
        node_on_visits = [[None for n in range(nb_customers+2)] for k in range(nb_trucks)]
        for k in range(nb_trucks):
            node_on_visits[k][0] = 0
            node_on_visits[k][nb_customers+1] = 0

        # During the route, the actual node visited is 1 + customers_sequences[k][n]:
        #  - When n < count(customers_sequences[k]) (when a customer is on the nth 
        #      position), node_on_visits[k][n+1] = cutomerId. (actual customers ids are 
        #      in [1..nb_customers], 0 being the warehouse)
        #  - When n >= count(customers_sequences[k]) (when the last customer has already 
        #      been visited), node_on_visits[k][n+1] = 0 (the warehouse id)
        for k in range(nb_trucks):
            for n in range(1, nb_customers+1):
                node_on_visits[k][n] = 1 + customers_sequences[k][n-1]

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

        # Create demands as an array to be able to access it with an "at" operator
        demands_array = model.array(demands)

        # The quantity needed in each route must not exceed the truck capacity
        route_quantity = [model.sum(demands_array[node_on_visits[k][n]] for n in range(1, nb_customers+1)) for k in range(nb_trucks)]
        for k in range(nb_trucks):
            model.constraint(route_quantity[k] <= truck_capacity)

        # Create distance as an array to be able to acces it with an "at" operator
        distance_array = model.array()
        for n in range(nb_customers+1):
            distance_array.add_operand(model.array(distance_matrix[n]))

        # Distance traveled by each truck
        route_distances = [None]*nb_trucks
        for k in range(nb_trucks):
            route_distances[k] = model.sum();
            for n in range(nb_customers+1):
                # Distance traveled from node n to node n+1 by truck k
                distance_node = model.at(distance_array, node_on_visits[k][n], node_on_visits[k][n+1])
                route_distances[k].add_operand(distance_node)

        # Total distance travelled
        total_distance = model.sum(route_distances)

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

        model.close()

        #
        # Parameterizes the solver
        #
        ls.create_phase().time_limit = int(str_time_limit)

        ls.solve()

        #
        # Writes the solution in a file with the following format :
        #  - number of trucks used and total distance
        #  - for each truck the nodes visited (omitting the start/end at the depot)
        #
        with open(sol_file, 'w') as f:
            f.write("%d %d\n" % (nb_trucks_used.value, total_distance.value))
            for k in range(nb_trucks):
                if(trucks_used[k].value != 1): continue
                # Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in [2..nbCustomers+1]
                # 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 "Augerat" format.
def read_input_cvrp(filename):
    file_it = iter(read_elem(sys.argv[1]))
    
    nb_nodes = 0
    while(1):
        token = file_it.next()
        if token == "DIMENSION":
            file_it.next() # Removes the ":"
            nb_nodes = int(file_it.next())
            nb_customers = nb_nodes - 1
        elif token == "CAPACITY":
            file_it.next() # Removes the ":"
            truck_capacity = int(file_it.next())
        elif token == "EDGE_WEIGHT_TYPE":
            file_it.next() # Removes the ":"
            token = file_it.next()
            if token != "EUC_2D":
                print ("Edge Weight Type " + token + " is not supported (only EUD_2D)")
                sys.exit(1)
        elif token == "NODE_COORD_SECTION":
            break;

    nodes_x = [None]*nb_nodes
    nodes_y = [None]*nb_nodes
    for n in range(nb_nodes):
        node_id = int(file_it.next())
        if node_id != n+1:
            print ("Unexpected index")
            sys.exit(1)
        nodes_x[n] = int(file_it.next())
        nodes_y[n] = int(file_it.next())

    # Compute distance matrix
    distance_matrix = compute_distance_matrix(nodes_x, nodes_y)

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

    demands = [None]*nb_nodes
    for n in range(nb_nodes):
        node_id = int(file_it.next())
        if node_id != n+1:
            print ("Unexpected index")
            sys.exit(1)
        demands[n] = int(file_it.next())

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

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

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

    if demands[0] != 0:
        print ("Warehouse demand is supposed to be 0")
        sys.exit(1)

    return (nb_customers, truck_capacity, distance_matrix, demands)


# Computes the distance matrix
def compute_distance_matrix(nodes_x, nodes_y):
    nb_nodes = len(nodes_x)
    distance_matrix = [[None for i in range(nb_nodes)] for j in range(nb_nodes)]
    for i in range(nb_nodes):
        distance_matrix[i][i] = 0
        for j in range(nb_nodes):
            dist = compute_dist(nodes_x[i], nodes_x[j], nodes_y[i], nodes_y[j])
            distance_matrix[i][j] = dist
            distance_matrix[j][i] = dist
    return distance_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))


def get_nb_trucks(filename):
    begin = filename.rfind("-k")
    if begin != -1:
        begin += 2
        end = filename.find(".", begin)
        return int(filename[begin:end])
    print ("Error: nb_trucks could not be read from the file name. Enter it from the command line")
    sys.exit(1)


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print ("Usage: python cvrp.py input_file [output_file] [time_limit] [nb_trucks]")
        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";
    str_nb_trucks = sys.argv[4] if len(sys.argv) > 4 else "0";

    main(instance_file, str_time_limit, sol_file, str_nb_trucks)
Compilation / Execution (Windows)
cl /EHsc cvrp.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver.dll.lib
cvrp instances\A-n32-k5.vrp
Compilation / Execution (Linux)
g++ cvrp.cpp -I/opt/localsolver_XXX/include -llocalsolver -lpthread -o cvrp
./cvrp instances/A-n32-k5.vrp
/********** cvrp.cpp **********/

#include <iostream>
#include <fstream>
#include <vector>
#include <cstring>
#include <cmath>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

class Cvrp {
public:
    // LocalSolver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each node
    vector<lsint> demands;

    // Distance matrix
    vector<vector<lsint> > distanceMatrix;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    vector<LSExpression> customersSequences;
    
    // Are the trucks actually used
    vector<LSExpression> trucksUsed;

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

    // Distance travelled by all the trucks
    LSExpression totalDistance;    
    
    /* Constructor */
    Cvrp(const char* strNbTrucks) {
        nbTrucks = atoi(strNbTrucks);
    }

    /* Reads instance data. */
    void readInstance(const string& fileName) {
        readInputCvrp(fileName);

        // The number of trucks is usually given in the name of the file
        // nbTrucks can also be given in command line
        if (nbTrucks == 0) nbTrucks = getNbTrucks(fileName);
    }

    void solve(int limit) {
        try {
            /* Declares the optimization model. */
            LSModel model = localsolver.getModel();

            // Sequence of customers visited by each truck.
            // When n >= count(customersSequences[k]) (after the last visited
            // customer), customersSequences[k][n] = -1
            customersSequences.resize(nbTrucks);
            for (int k = 0; k < nbTrucks; k++) {
                customersSequences[k] = model.listVar(nbCustomers);
            }

            // All customers must be visited by  the trucks
            model.constraint(model.partition(customersSequences.begin(), customersSequences.end()));

            // Each truck starts (n=0) and ends (n=nbCustomers+1) at the depot (index 0)
            vector<vector<LSExpression> > nodeOnVisits(nbTrucks);
            for (int k = 0; k < nbTrucks; k++) {
                nodeOnVisits[k].resize(nbCustomers+2);
                nodeOnVisits[k][0] = model.createConstant((lsint) 0);
                nodeOnVisits[k][nbCustomers+1] = model.createConstant((lsint) 0);
            }

            // During the route, the actual node visited is 1 + customersSequences[k][n]:
            //  - When n < count(customersSequences[k]) (when a customer is on the nth 
            //      position), nodeOnVisits[k][n+1] = cutomerId. (actual customers ids are 
            //      in [1..nbCustomers], 0 being the warehouse)
            //  - When n >= count(customersSequences[k]) (when the last customer has already 
            //      been visited), nodeOnVisits[k][n+1] = 0 (the warehouse id)
            for (int k = 0; k < nbTrucks; k++) {
                for (int n = 1; n < nbCustomers + 1; n++) {
                    nodeOnVisits[k][n] = 1 + customersSequences[k][n-1];
                }
            }

            // A truck is used if it visits at least one customer
            trucksUsed.resize(nbTrucks);
            for (int k = 0; k < nbTrucks; k++) {
                trucksUsed[k] = model.count(customersSequences[k]) > 0;
            }
            nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());

            // Create demands as an array to be able to access it with an "at" operator
            LSExpression demandsArray = model.array(demands.begin(), demands.end());

            // The quantity needed in each route must not exceed the truck capacity
            for (int k = 0; k < nbTrucks; k++) {
                LSExpression routeQuantity = model.sum();
                for (int n = 1; n < nbCustomers + 1; n++) {
                    routeQuantity += demandsArray[nodeOnVisits[k][n]];
                }
                model.constraint(routeQuantity <= truckCapacity);
            }

            // Create distance as an array to be able to acces it with an "at" operator
            LSExpression distanceArray = model.array();
            for (int n = 0; n < nbCustomers + 1; n++) {
                distanceArray.addOperand(model.array(distanceMatrix[n].begin(), distanceMatrix[n].end()));
            }

            // Distance traveled by each truck
            vector<LSExpression> routeDistances(nbTrucks);
            for (int k = 0; k < nbTrucks; k++) {
                routeDistances[k] = model.sum();
                for (int n = 0; n < nbCustomers + 1; n++) {
                    // Distance traveled from node n to node n+1 by truck k
                    routeDistances[k] += model.at(distanceArray, nodeOnVisits[k][n], nodeOnVisits[k][n+1]);
                }
            }

            // Total distance travelled
            totalDistance = model.sum(routeDistances.begin(), routeDistances.end());

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

            model.close();

            /* Parameterizes the solver. */
            LSPhase phase = localsolver.createPhase();
            phase.setTimeLimit(limit);
            localsolver.solve();
            
        } catch (const LSException& e) {
            cout << "LSException:" << e.getMessage() << endl;
            exit(1);
        }
    }

    /* Writes the solution in a file with the following format :
     *  - number of trucks used and total distance
     *  - for each truck the nodes visited (omitting the start/end at the depot)*/
    void writeSolution(const string& fileName) {
        ofstream outfile(fileName.c_str());
        if (!outfile.is_open()) {
            cerr << "File " << fileName << " cannot be opened." << endl;
            exit(1);
        }
        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-1]. +2 is to put it back in [2..nbCustomers+1]
            // as in the data files (1 being the depot)
            LSCollection customersCollection = customersSequences[k].getCollectionValue();
            for (lsint i = 0; i < customersCollection.count(); i++) {
                outfile << customersCollection[i] + 2 << " ";
            }
            outfile << endl;
        }
        outfile.close();
    }

private:
    /* The input files follow the "Augerat" format.*/
    void readInputCvrp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            cerr << "File " << fileName << " cannot be opened." << endl;
            exit(1);
        }
        
        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, "CAPACITY") == 0) {
                pch = strtok(NULL, " :");
                truckCapacity = atoi(pch);
            } else if (strcmp(pch, "EDGE_WEIGHT_TYPE") == 0) {
                pch = strtok(NULL, " :");
                if (strcmp(pch, "EUC_2D") != 0) {
                    cerr << "Edge Weight Type " << pch << " is not supported (only EUC_2D)" << endl;
                    exit(1);
                }
            } else if (strcmp(pch, "NODE_COORD_SECTION") == 0) {
                break;
            }
        }

        vector<int> nodesX(nbNodes);
        vector<int> nodesY(nbNodes);
        for (int n = 0; n < nbNodes; n++) {
            int id;
            infile >> id;
            if (id != n+1) {
                cerr << "Unexpected index" << endl;
                exit(1);
            }
            infile >> nodesX[n];
            infile >> nodesY[n];
        }

        // Compute distance matrix
        computeDistanceMatrix(nodesX, nodesY);

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " :");
        if (strcmp(pch, "DEMAND_SECTION") != 0) {
            cerr << "Expected keyword DEMAND_SECTION" << endl;
            exit(1);
        }

        demands.resize(nbNodes);
        for (int n = 0; n < nbNodes; n++) {
            int id;
            infile >> id;
            if (id != n+1) {
                cerr << "Unexpected index" << endl;
                exit(1);
            }
            infile >> demands[n];
        }

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " :");
        if (strcmp(pch, "DEPOT_SECTION") != 0) {
            cerr << "Expected keyword DEPOT_SECTION" << endl;
            exit(1);
        }

        int warehouseId;
        infile >> warehouseId;
        if (warehouseId != 1) {
            cerr << "Warehouse id is supposed to be 1" << endl;
            exit(1);
        }

        int endOfDepotSection;
        infile >> endOfDepotSection;
        if (endOfDepotSection != -1) {
            cerr << "Expecting only one warehouse, more than one found" << endl;
            exit(1);
        }

        if (demands[0] != 0) {
            cerr << "Warehouse demand is supposed to be 0" << endl;
            exit(1);
        }

        infile.close();
    }

    /* Computes the distance matrix */
    void computeDistanceMatrix(const vector<int>& nodesX, const vector<int>& nodesY) {
        distanceMatrix.resize(nbCustomers+1);
        for (int i = 0; i < nbCustomers+1; i++) {
            distanceMatrix[i].resize(nbCustomers+1);
        }
        for (int i = 0; i < nbCustomers+1; i++) {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers + 1; j++) {
                lsint distance = computeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrix[i][j] = distance;
                distanceMatrix[j][i] = distance;
            }
        }
    }

    lsint 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 getNbTrucks(const string& fileName) {
        size_t pos = fileName.rfind("-k");
        if (pos != string::npos) {
            string nbTrucksStr = fileName.substr(pos + 2);
            pos = nbTrucksStr.find(".");
            if (pos != string::npos) {
                return atoi(nbTrucksStr.substr(0, pos).c_str());
            }
        }
        cerr << "Error: nbTrucks could not be read from the file name. Enter it from the command line" << endl;
        exit(1);
        return -1;
    }
};

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

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

    Cvrp model(strNbTrucks);
    model.readInstance(instanceFile);
    model.solve(atoi(strTimeLimit));
    if(solFile != NULL)
        model.writeSolution(solFile);
    return 0;
}
Compilation / Execution (Windows)
javac Cvrp.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. Cvrp instances\A-n32-k5.vrp
Compilation/Execution (Linux)
javac Cvrp.java -cp /opt/localsolver_XXX/bin/localsolver.jar
java -cp /opt/localsolver_XXX/bin/localsolver.jar:. Cvrp instances/A-n32-k5.vrp
/********** Cvrp.java **********/

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

public class Cvrp {
    // Solver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each node
    long[] demands;

    // Distance matrix
    long[][] distanceMatrix;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    LSExpression[] customersSequences;
    
    // Are the trucks actually used
    LSExpression[] trucksUsed;

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

    // Distance travelled by all the trucks
    LSExpression totalDistance;

    Cvrp(String strNbTrucks) {
        localsolver = new LocalSolver();
        nbTrucks = Integer.parseInt(strNbTrucks);
    }

    /* Reads instance data. */
    void readInstance(String fileName) {
        readInputCvrp(fileName);

        // The number of trucks is usually given in the name of the file
        // nbTrucks can also be given in command line
        if (nbTrucks == 0) nbTrucks = getNbTrucks(fileName);
    }

    void solve(int limit) {
        try {
            /* Declares the optimization model. */
            LSModel model = localsolver.getModel();

            // Sequence of customers visited by each truck.
            // When n >= count(customersSequences[k]) (after the last visited
            // customer), customersSequences[k][n] = -1
            customersSequences = new LSExpression[nbTrucks];
            for (int k = 0; k < nbTrucks; k++) {
                customersSequences[k] = model.listVar(nbCustomers);
            }

            // All customers must be visited by  the trucks
            model.constraint(model.partition(customersSequences));

            // Each truck starts (n=0) and ends (n=nbCustomers+1) at the depot (index 0)
            LSExpression[][] nodeOnVisits = new LSExpression[nbTrucks][nbCustomers+2];
            for (int k = 0; k < nbTrucks; k++) {
                nodeOnVisits[k][0] = model.createConstant(0);
                nodeOnVisits[k][nbCustomers+1] = model.createConstant(0);
            }

            // During the route, the actual node visited is 1 + customersSequences[k][n]:
            //  - When n < count(customersSequences[k]) (when a customer is on the nth 
            //      position), nodeOnVisits[k][n+1] = cutomerId. (actual customers ids are 
            //      in [1..nbCustomers], 0 being the warehouse)
            //  - When n >= count(customersSequences[k]) (when the last customer has already 
            //      been visited), nodeOnVisits[k][n+1] = 0 (the warehouse id)
            for (int k = 0; k < nbTrucks; k++) {
                for (int n = 1; n < nbCustomers + 1; n++) {
                    LSExpression customerOnVisit = model.at(customersSequences[k], n-1);
                    nodeOnVisits[k][n] = model.sum(1, customerOnVisit);
                }
            }

            // A truck is used if it visits at least one customer
            trucksUsed = new LSExpression[nbTrucks];
            for (int k = 0; k < nbTrucks; k++) {
                trucksUsed[k] = model.gt(model.count(customersSequences[k]), 0);
            }
            nbTrucksUsed = model.sum(trucksUsed);

            // Create demands as an array to be able to access it with an "at" operator
            LSExpression demandsArray = model.array(demands);

            // The quantity needed in each route must not exceed the truck capacity
            for (int k = 0; k < nbTrucks; k++) {
                LSExpression routeQuantity = model.sum();
                for (int n = 1; n < nbCustomers + 1; n++) {
                    routeQuantity.addOperand(model.at(demandsArray, nodeOnVisits[k][n]));
                }
                model.constraint(model.leq(routeQuantity, truckCapacity));
            }

            // Create distance as an array to be able to acces it with an "at" operator
            LSExpression distanceArray = model.array();
            for (int n = 0; n < nbCustomers + 1; n++) {
                distanceArray.addOperand(model.array(distanceMatrix[n]));
            }

            // Distance traveled by each truck
            LSExpression[] routeDistances = new LSExpression[nbTrucks];
            for (int k = 0; k < nbTrucks; k++) {
                routeDistances[k] = model.sum();
                for (int n = 0; n < nbCustomers + 1; n++) {
                    // Distance traveled from node n to node n+1 by truck k
                    routeDistances[k].addOperand(model.at(distanceArray, nodeOnVisits[k][n], nodeOnVisits[k][n+1]));
                }
            }

            // Total distance travelled
            totalDistance = model.sum(routeDistances);

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

            model.close();

            /* Parameterizes the solver. */
            LSPhase phase = localsolver.createPhase();
            phase.setTimeLimit(limit);

            localsolver.solve();

        } catch (LSException e) {
            System.out.println("LSException:" + e.getMessage());
            System.exit(1);
        }
    }

    /* Writes the solution in a file with the following format :
     *  - number of trucks used and total distance
     *  - for each truck the nodes visited (omitting the start/end at the depot)*/
    void writeSolution(String fileName) {
        try {
            BufferedWriter output = new BufferedWriter(new FileWriter(fileName));
            output.write(nbTrucksUsed.getValue() + " " + totalDistance.getValue() + "\n");
            for (int k = 0; k < nbTrucks; k++) {
                if (trucksUsed[k].getValue() != 1) continue;
                // Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in [2..nbCustomers+1]
                // as in the data files (1 being the depot)
                LSCollection customersCollection = customersSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); i++) {
                    output.write((customersCollection.get(i) + 2) + " ");
                }
                output.write("\n");
            }
            output.close();
        } catch (IOException ex) {
            ex.printStackTrace();
        }
    }
    
     public static void main(String[] args) {
        if (args.length < 1) {
            System.out.println("Usage: java Cvrp inputFile [outputFile] [timeLimit] [nbTrucks]");
            System.exit(1);
        }

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

        Cvrp model = new Cvrp(strNbTrucks);
        model.readInstance(instanceFile);
        model.solve(Integer.parseInt(strTimeLimit));
        if(outputFile != null) {
            model.writeSolution(outputFile);
        }
    }

    /* The input files follow the "Augerat" format.*/
    private void readInputCvrp(String fileName) { 
        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("CAPACITY")) {
                    truckCapacity = Integer.parseInt(splitted[1].trim());
                }
                else if (splitted[0].contains("EDGE_WEIGHT_TYPE")) {
                    if (splitted[1].trim().compareTo("EUC_2D") != 0) {
                        System.out.println("Edge Weight Type " + splitted[1] + " is not supported (only EUC_2D)");
                        System.exit(1);
                    }
                }
                else if (splitted[0].contains("NODE_COORD_SECTION")) {
                    break;
                }
            }

            int[] nodesX = new int[nbNodes];
            int[] nodesY = new int[nbNodes];
            for (int n = 0; n < nbNodes; n++) {
                int id = input.nextInt();
                if (id != n + 1) {
                    System.out.println("Unexpected index");
                    System.exit(1);
                }
                nodesX[n] = input.nextInt();
                nodesY[n] = input.nextInt();
            }

            computeDistanceMatrix(nodesX, nodesY);

            splitted = input.nextLine().split(":"); //End the last line
            splitted = input.nextLine().split(":");
            if (!splitted[0].contains("DEMAND_SECTION")) {
                System.out.println("Expected keyword DEMAND_SECTION");
                System.exit(1);
            }

            demands = new long[nbNodes];
            for (int n = 0; n < nbNodes; n++) {
                int id = input.nextInt();
                if (id != n + 1) {
                    System.out.println("Unexpected index");
                    System.exit(1);
                }
                demands[n] = input.nextInt();
            }

            splitted = input.nextLine().split(":"); //End the last line
            splitted = input.nextLine().split(":");
            if (!splitted[0].contains("DEPOT_SECTION")) {
                System.out.println("Expected keyword DEPOT_SECTION");
                System.exit(1);
            }

            int warehouseId = input.nextInt();
            if (warehouseId != 1) {
                System.out.println("Warehouse id is supposed to be 1");
                System.exit(1);
            }

            int endOfDepotSection = input.nextInt();
            if (endOfDepotSection != -1) {
                System.out.println("Expecting only one warehouse, more than one found");
                System.exit(1);
            }

            if (demands[0] != 0) {
                System.out.println("Warehouse demand is supposed to be 0");
                System.exit(1);
            }

            input.close();
            
        } catch (IOException ex) {
            ex.printStackTrace();
        }
    }

    /* Computes the distance matrix */
    private void computeDistanceMatrix(int[] nodesX, int[] nodesY) {
        distanceMatrix = new long[nbCustomers+1][nbCustomers+1];
        for (int i = 0; i < nbCustomers + 1; i++) {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers + 1; j++) {
                long dist = computeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrix[i][j] = dist;
                distanceMatrix[j][i] = dist;
            }
        }
    }

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

    private int getNbTrucks(String fileName) {
        int begin = fileName.lastIndexOf("-k");
        if (begin != -1) {
            begin += 2;
            int end = fileName.indexOf(".", begin);
            return Integer.parseInt(fileName.substring(begin, end));
        }
        System.out.println("Error: nbTrucks could not be read from the file name. Enter it from the command line");
        System.exit(1);
        return -1;
    }
}
Compilation/Execution (Windows)
copy %LS_HOME%\bin\*net.dll .
csc Cvrp.cs /reference:localsolvernet.dll
Cvrp instances\A-n32-k5.vrp
/********** Cvrp.cs **********/

using System;
using System.IO;
using localsolver;


public class Cvrp : IDisposable
{
    // Solver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each node
    long[] demands;

    // Distance matrix
    long[][] distanceMatrix;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    LSExpression[] customersSequences;
    
    // Are the trucks actually used
    LSExpression[] trucksUsed;

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

    // Distance travelled by all the trucks
    LSExpression totalDistance;

    public Cvrp (string strNbTrucks)
    {
        localsolver = new LocalSolver();
        nbTrucks = int.Parse(strNbTrucks);
    }

    /* Reads instance data. */
    void ReadInstance(string fileName)
    {
        readInputCvrp(fileName);

        // The number of trucks is usually given in the name of the file
        // nbTrucks can also be given in command line
        if (nbTrucks == 0) nbTrucks = getNbTrucks(fileName);
    }

    public void Dispose ()
    {
        localsolver.Dispose();
    }

    void Solve(int limit)
    {
        /* Declares the optimization model. */
        LSModel model = localsolver.GetModel();

        // Sequence of customers visited by each truck.
        // When n >= count(customersSequences[k]) (after the last visited
        // customer), customersSequences[k][n] = -1
        customersSequences = new LSExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; k++)
        {
            customersSequences[k] = model.List(nbCustomers);
        }

        // All customers must be visited by  the trucks
        model.Constraint(model.Partition(customersSequences));

        // Each truck starts (n=0) and ends (n=nbCustomers+1) at the depot (index 0)
        LSExpression[,] nodeOnVisits = new LSExpression[nbTrucks,nbCustomers+2];
        for (int k = 0; k < nbTrucks; k++)
        {
            nodeOnVisits[k,0] = model.CreateConstant(0);
            nodeOnVisits[k,nbCustomers+1] = model.CreateConstant(0);
        }

        // During the route, the actual node visited is 1 + customersSequences[k][n]:
        //  - When n < count(customersSequences[k]) (when a customer is on the nth 
        //      position), nodeOnVisits[k][n+1] = cutomerId. (actual customers ids are 
        //      in [1..nbCustomers], 0 being the warehouse)
        //  - When n >= count(customersSequences[k]) (when the last customer has already 
        //      been visited), nodeOnVisits[k][n+1] = 0 (the warehouse id)
        for (int k = 0; k < nbTrucks; k++)
        {
            for (int n = 1; n < nbCustomers + 1; n++)
            {
                nodeOnVisits[k,n] = 1 + customersSequences[k][n-1];
            }
        }

        // A truck is used if it visits at least one customer
        trucksUsed = new LSExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; k++)
        {
            trucksUsed[k] = model.Count(customersSequences[k]) > 0;
        }
        nbTrucksUsed = model.Sum(trucksUsed);

        // Create demands as an array to be able to access it with an "at" operator
        LSExpression demandsArray = model.Array(demands);

        // The quantity needed in each route must not exceed the truck capacity
        for (int k = 0; k < nbTrucks; k++)
        {
            LSExpression routeQuantity = model.Sum();
            for (int n = 1; n < nbCustomers + 1; n++)
            {
                routeQuantity.AddOperand(demandsArray[nodeOnVisits[k,n]]);
            }
            model.Constraint(routeQuantity <= truckCapacity);
        }

        // Create distance as an array to be able to acces it with an "at" operator
        LSExpression distanceArray = model.Array();
        for (int n = 0; n < nbCustomers + 1; n++)
        {
            distanceArray.AddOperand(model.Array(distanceMatrix[n]));
        }

        // Distance traveled by each truck
        LSExpression[] routeDistances = new LSExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; k++)
        {
            routeDistances[k] = model.Sum();
            for (int n = 0; n < nbCustomers + 1; n++)
            {
                // Distance traveled from node n to node n+1 by truck k
                // distanceArray[a, b] is a shortcut for accessing the multi-dimensional array
                // distanceArray with an at operator. Same as model.At(distanceArray, a, b)
                routeDistances[k].AddOperand(distanceArray[nodeOnVisits[k,n], nodeOnVisits[k,n+1]]);
            }
        }

        // Total distance travelled
        totalDistance = model.Sum(routeDistances);

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

        model.Close();
        
        /* Parameterizes the solver. */
        LSPhase phase = localsolver.CreatePhase();
        phase.SetTimeLimit(limit);

        localsolver.Solve();
    }

    /* Writes the solution in a file with the following format :
     *  - number of trucks used and total distance
     *  - for each truck the nodes 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-1]. +2 is to put it back in [2..nbCustomers+1]
                // as in the data files (1 being the depot)
                LSCollection 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: Cvrp inputFile [solFile] [timeLimit] [nbTrucks]");
            System.Environment.Exit (1);
        }
        String instanceFile = args [0];
        String outputFile = args.Length > 1 ? args [1] : null;
        String strTimeLimit = args.Length > 2 ? args [2] : "20";
        String strNbTrucks = args.Length > 3 ? args [3] : "0";

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

    /* The input files follow the "Augerat" format.*/
    private void readInputCvrp(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("CAPACITY"))
                {
                    truckCapacity = int.Parse(splitted[1]);
                }
                else if (splitted[0].Contains("EDGE_WEIGHT_TYPE"))
                {
                    if (!splitted[1].Trim().Equals("EUC_2D"))
                    {
                        Console.WriteLine("Edge Weight Type " + splitted[1] + " is not supported (only EUC_2D)");
                        System.Environment.Exit(1);
                    }
                }
                else if (splitted[0].Contains("NODE_COORD_SECTION"))
                {
                    break;
                }
            }
            int[] nodesX = new int[nbNodes];
            int[] nodesY = new int[nbNodes];
            for (int n = 0; n < nbNodes; n++)
            {
                splitted = input.ReadLine().Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n+1)
                {
                    Console.WriteLine("Unexpected index");
                    System.Environment.Exit(1);
                }
                nodesX[n] = int.Parse(splitted[1]);
                nodesY[n] = int.Parse(splitted[2]);
            }

            computeDistanceMatrix(nodesX, nodesY);

            splitted = input.ReadLine().Split(':');
            if (!splitted[0].Contains("DEMAND_SECTION"))
            {
                Console.WriteLine("Expected keyword DEMAND_SECTION");
                System.Environment.Exit(1);
            }

            demands = new long[nbNodes];
            for (int n = 0; n < nbNodes; n++)
            {
                splitted = input.ReadLine().Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n+1)
                {
                    Console.WriteLine("Unexpected index");
                    System.Environment.Exit(1);
                }
                demands[n] = int.Parse(splitted[1]);
            }

            splitted = input.ReadLine().Split(':');
            if (!splitted[0].Contains("DEPOT_SECTION"))
            {
                Console.WriteLine("Expected keyword DEPOT_SECTION");
                System.Environment.Exit(1);
            }

            int warehouseId = int.Parse(input.ReadLine());
            if (warehouseId != 1)
            {
                Console.WriteLine("Warehouse id is supposed to be 1");
                System.Environment.Exit(1);
            }

            int endOfDepotSection = int.Parse(input.ReadLine());
            if (endOfDepotSection != -1)
            {
                Console.WriteLine("Expecting only one warehouse, more than one found");
                System.Environment.Exit(1);
            }

            if (demands[0] != 0)
            {
                Console.WriteLine("Warehouse demand is supposed to be 0");
                System.Environment.Exit(1);
            }
        }
    }

    /* Computes the distance matrix */
    private void computeDistanceMatrix(int[] nodesX, int[] nodesY)
    {
        distanceMatrix = new long[nbCustomers+1][];
        for (int i = 0; i < nbCustomers + 1; i++)
        {
            distanceMatrix[i] = new long[nbCustomers+1];
        }
        for (int i = 0; i < nbCustomers + 1; i++)
        {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers + 1; j++)
            {
                long dist = computeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrix[i][j] = dist;
                distanceMatrix[j][i] = dist;
            }
        }
    }

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

    private int getNbTrucks(string fileName)
    {
        string[] splitted = fileName.Split(new Char[] {'-', 'k'}, StringSplitOptions.RemoveEmptyEntries);
        if (splitted.Length >= 2)
        {
            String toSplit = splitted[splitted.Length - 1];
            splitted = toSplit.Split(new Char[] {'.'}, StringSplitOptions.RemoveEmptyEntries);
            return int.Parse(splitted[0]);
        }
        Console.WriteLine("Error: nbTrucks could not be read from the file name. Enter it from the command line");
        System.Environment.Exit(1);
        return -1;
    }
}

From CVRP to CVRPTWΒΆ

The Capacitated Vehicule Routing Problem with Time Windows (CVRPTW) is a variant of the CVRP we just modelled. In addition to the vehicule capacity constraint, each customer has a service time and a time window during which it can be serviced.

To take the time windows into account in our CVRP model, we have to track the arrival time of the trucks on each visit. If we do not take into account the start of the time windows, this arrival time is the arrival time on the last visit, plus the service time there, plus the travel time from the last visit to the current. However, when a truck arrives at a customer, it has to wait the beginning of the time window in case it is too early. Let’s assume we have three LocalSolver arrays: startTW, endTW and serviceTime, corresponding respectively to the start and the end of the time window and the service time of each node. The arrival time would be defined as follows:

arrivalTimes[k in 1..nbTrucks][0] <- startTW[0]; // Start at the opening of the depot
arrivalTimes[k in 1..nbTrucks][n in 1..nbCustomers+1]
    <- max(startTW[nodeOnVisits[k][n]],
           arrivalTimes[k][n-1] + serviceTime[nodeOnVisits[k][n-1]] + distanceNodes[k][n-1]);

The violation of the end of a customer time window is more a first-priority objective than a real constraint (cf. this article for more information). Indeed, a customer would often wait after the end of its time window if the delivery still has not occured. Only afterwards would he eventually point out to the company that his time-window has not been respected and he now needs compensation. On the contrary, the respect of the horizon is a structural constraint: trucks must respect the closing hours of the depot (its end of time window) as it can be regulated by laws or by company rules. It can be done as follows:

violations[k in 1..nbTrucks][n in 1..nbCustomers]
    <- max(arrivalTimes[k][n] - endTW[nodeOnVisits[k][n]], 0);
// Constraint on the horizon
for[k in 1..nbTrucks]
    constraint arrivalTimes[k][nbCustomers+1] <= endTW[0];

Warning

Along each tour, after the last visit n*, the truck goes to the depot on visit n*+1 and then does virtual visits from the depot to the depot for each n in [n*+1..nbCustomers+1] (because nodeOnVisits[k][n] is 0). If the service time of the depot is not 0, the arrival time on each of these virtual visit will increment from this service time. As a consequence, the last arrival time to the depot arrivalTimes[k][nbCustomers+1] will be greater than the actual one arrivalTimes[k][n*+1]. As we are in a context where trucks start and end at the depot with no possible refill in between, a depot service time only means a loading time of the trucks at the beginning. Provided that when we read the instance we saved the loading time into loadingTime and that we took care of setting serviceTime[0] to 0, it is possible to take it into account simply by changing the arrival time on visit 0 as follows:

arrivalTimes[k in 1..nbTrucks][0] <- startTW[0] + loadingTime;

The model will now have 3 objectives to minimize: the sum of violations, the number of trucks and the total distance. It could be hard for LocalSolver to jump from a solution with no violations to another solution with a truck less because it would also need to have no violations (because of the lexicographical order). These two solutions can be really far away from each other or even have no path linking them. In order to efficiently solve the CVRPTW problem, a good practice is to solve the problem step by step:

  • Construct a model with only the violations objective and a fixed number of trucks. LocalSolver will return optimal when the violation objective reaches 0 (its natural lower bound).

  • Solve this model following this pseudo-code:

    nbTrucks = n
    lastSolution = None
    while(true)
        solve cvrptw nbTrucks-1 timeLimit
        if solution is optimal
            nbTrucks--
            lastSolution = solution
        else
            break
    
  • Retrieve lastSolution and nbTrucks

  • Add the distance objective to the model (it now has 2 objectives: first minimize the violations, then minimize the distance)

  • Solve this model with nbTrucks and the initial solution set to lastSolution

It can be done whether through the APIs or also in LSP language, with a bit of bash scripting on top of it. It is a very good practice to apply to this kind of problems. It also proved to be really efficient with binpacking problems: solve iteratively the packing problem with one bin less each time.