Vehicle routing problem

Principles learned

  • Add multiple list decision variables
  • Add a partition constraint
  • Use lambda expression to compute a sum with a variable number of terms
  • 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 vehicles with uniform capacity must service customers with known demand for a single commodity. The vehicles start and end their routes at a common depot. Each customer can only be served by one vehicle. 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 truck 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 a 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.

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. The total quantity delivered by each truck is computed with a function to apply the sum operator over all customers visited by a truck. Note that the number of terms in this sum varies automatically with the size of the list. This quantity is constrained to be lower than the truck capacity.

For each truck, the distance travelled from the visit n-1 to visit n is accessed with an “at” operator to the multi-dimensional array distanceMatrix, with the first index. Here again we use a function to sum distances along each route.

The two 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.

If you are interested in the classical variant where time-windows are defined for each customer, you can now study our CVRPTW model.

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.
    customersSequences[k in 1..nbTrucks] <- list(nbCustomers);

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

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

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

        // The quantity needed in each route must not exceed the truck capacity
        routeQuantity <- sum(0..c-1, i => demands[sequence[i]]);
        constraint routeQuantity <= truckCapacity;

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

    nbTrucksUsed <- sum[k in 1..nbTrucks](trucksUsed[k]);

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

    //nodeX and nodeY are indexed by original data indices (1 for depot)
    for[n in 1..nbNodes] {
        if (n != inFile.readInt()) throw "Unexpected index";
        nodesX[n] = round(inFile.readDouble());
        nodesY[n] = 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.
        local demand = inFile.readInt();
        if (n == 1) {
            if (demand != 0) error("expected demand for depot is 0");
        } else {
            demands[n-2] = demand; // demands is indexed by customers
        }
    }

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

/* Compute the distance between each node */
function computeDistanceMatrix() {
    for[i in 0..nbCustomers-1]
    {
        distanceMatrix[i][i] = 0;
        for[j in i+1..nbCustomers-1]
        {        
            // +2 because computeDist expected original data indices,
            // with customers in 2..nbNodes (1 being the depot)
            local localDistance = computeDist(i+2, j+2);
            distanceMatrix[j][i] = localDistance;
            distanceMatrix[i][j] = localDistance;
        }
    }

    for[i in 0..nbCustomers-1] {
        local localDistance = computeDist(1, i+2);
        distanceWarehouse[i] = 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, distance_warehouses, 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.
        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))
        
         # Create demands as an array to be able to access it with an "at" operator
        demands_array = model.array(demands)
        
         # 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):
            distance_array.add_operand(model.array(distance_matrix[n]))

        distance_warehouse_array = model.array(distance_warehouses)

        route_distances = [None for n in range(nb_trucks)]

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

        for k in range(nb_trucks):
            sequence = customers_sequences[k]
            c = model.count(sequence)
                        
            # Quantity in each truck
            demand_selector = model.function(lambda i: model.at(demands_array, sequence[i]))
            route_quantity = model.sum(model.range(0, c), demand_selector) 
            model.constraint(route_quantity <= truck_capacity)
            
            # Distance traveled by each truck
            dist_selector = model.function(lambda i: model.at(distance_array, sequence[i-1], sequence[i]))
            route_distances[k] = model.sum(model.range(1,c), dist_selector) + \
                 model.iif(c > 0, model.at(distance_warehouse_array, sequence[0]) + model.at(distance_warehouse_array, sequence[c-1]),0)                       
       
        # 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)
        #
        if len(sys.argv) >= 3:
            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;

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

    # Compute distance matrix
    distance_matrix = compute_distance_matrix(customers_x, customers_y)
    distance_warehouses = compute_distance_warehouses(depot_x, depot_y, customers_x, customers_y)

    token = file_it.next()
    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(file_it.next())
        if node_id != n+1:
            print ("Unexpected index")
            sys.exit(1)
        if node_id == 1:
            if  int(file_it.next()) != 0:
                print ("Demand for depot should be 0")
                sys.exit(1)
        else:   
            # -2 because orginal customer indices are in 2..nbNodes
            demands[node_id-2] = 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)

    return (nb_customers, truck_capacity, distance_matrix, distance_warehouses, demands)


# Computes 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


# Computes the distances to warehouse
def compute_distance_warehouses(depot_x, depot_y, customers_x, customers_y):
    nb_customers = len(customers_x)
    distance_warehouses = [None for i in range(nb_customers)]
    for i in range(nb_customers):
        dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
        distance_warehouses[i] = dist
    return distance_warehouses



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;

    // Distance to depot
    vector<lsint> distanceWarehouses;

    // 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) {
            // Declares the optimization model.
            LSModel model = localsolver.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  the trucks
            model.constraint(model.partition(customersSequences.begin(), customersSequences.end()));

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

            // 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; n++) {
                distanceArray.addOperand(model.array(distanceMatrix[n].begin(), distanceMatrix[n].end()));
            }
            LSExpression distanceWarehousesArray = model.array(distanceWarehouses.begin(), distanceWarehouses.end());
            
            trucksUsed.resize(nbTrucks);
            vector<LSExpression> routeDistances(nbTrucks);
                       
            for (int k = 0; k < nbTrucks; k++) {
                LSExpression sequence = customersSequences[k];
                LSExpression c = model.count(sequence);

                // A truck is used if it visits at least one customer
                trucksUsed[k] = c > 0;
                                
                // The quantity needed in each route must not exceed the truck capacity
		        LSExpression demandSelector = model.createFunction([&](LSExpression i) { return demandsArray[sequence[i]]; });	
		        LSExpression routeQuantity = model.sum(model.range(0, c), demandSelector);
		        model.constraint(routeQuantity <= truckCapacity);

		        // Distance traveled by truck k
		        LSExpression distSelector = model.createFunction([&](LSExpression i) { return model.at(distanceArray, sequence[i - 1], sequence[i]); });	 
                routeDistances[k] = model.sum(model.range(1, c), distSelector) + 
                    model.iif(c > 0, distanceWarehousesArray[sequence[0]] + distanceWarehousesArray[sequence[c - 1]], 0);
            }
            
            // Total nb trucks used
            nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());

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

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

private:
    // The input files follow the "Augerat" format.
    void readInputCvrp(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, "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) {
                    throw std::runtime_error("Only Edge Weight Type EUC_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 orginal customer indices are in 2..nbNodes 
                infile >> customersX[n-2];
                infile >> customersY[n-2];
            }
        }

        // Compute distance matrix
        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");
        }

        demands.resize(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 orginal customer indices are in 2..nbNodes 
                demands[n-2] = demand;
            }                
        }

        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 warehouseId;
        infile >> warehouseId;
        if (warehouseId != 1) {
            throw std::runtime_error("Warehouse id is supposed to be 1");
        }

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

        infile.close();
    }

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

    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());
            }
        }
        throw std::runtime_error("Error: nbTrucks could not be read from the file name. Enter it from the command line");
        return -1;
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: cvrp inputFile [outputFile] [timeLimit] [nbTrucks]" << endl;
        return 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";

    try {
        Cvrp model(strNbTrucks);
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if(solFile != NULL) model.writeSolution(solFile);
        return 0;
    } catch (const exception& e){
        cerr << "Error occurred: " << e.what() << endl;
        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 (= number of nodes minus 1)
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each customer
    long[] demands;

    // Distance matrix between customers
    long[][] distanceMatrix;

    // Distances between customers and warehouse
    long[] distanceWarehouses;

    // Number of trucks
    int nbTrucks;

    // Decision variables
    LSExpression[] customersSequences;

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

    // Distance travelled by each truck
    LSExpression[] routeDistances;

    // 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()
    {
        if (localsolver != null)
            localsolver.Dispose();
    }

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

        trucksUsed = new LSExpression[nbTrucks];
        customersSequences = new LSExpression[nbTrucks];
        routeDistances = new LSExpression[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 the trucks
        model.Constraint(model.Partition(customersSequences));

        // Create demands and distances as arrays to be able to access it with an "at" operator
        LSExpression demandsArray = model.Array(demands);
        LSExpression distanceWarehouseArray = model.Array(distanceWarehouses);
        LSExpression distanceArray = model.Array(distanceMatrix);

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

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

            // The quantity needed in each route must not exceed the truck capacity
            LSExpression demandSelector = model.Function(i => demandsArray[sequence[i]]);
            LSExpression routeQuantity = model.Sum(model.Range(0, c), demandSelector);
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance travelled by truck k
            LSExpression distSelector = model.Function(i => distanceArray[sequence[i - 1], sequence[i]]);
            routeDistances[k] = model.Sum(model.Range(1, c), distSelector)
                                + model.If(c > 0, distanceWarehouseArray[sequence[0]] + distanceWarehouseArray[sequence[c - 1]], 0);
        }

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

        // Objective: minimize the number of trucks used, then minimize the distance traveled
        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]");
            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"))
                        throw new Exception("Edge Weight Type " + splitted[1] + " is not supported (only EUC_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 orginal 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");

            demands = 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("Warehouse demand is supposed to be 0");
                } else {
                    // -2 because orginal customer indices are in 2..nbNodes 
                    demands[n-2] = demand;
                }
            }

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

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

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

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

        for (int i = 0; i < nbCustomers; i++)
        {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; j++)
            {
                long dist = ComputeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distanceMatrix[i][j] = dist;
                distanceMatrix[j][i] = dist;
            }
        }

        distanceWarehouses = new long[nbCustomers];
        for(int i = 0; i < nbCustomers; ++i) {
            distanceWarehouses[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));
    }

    private int GetNbTrucks(string fileName)
    {
        string[] splitted = fileName.Split(new[] { '-', 'k' }, StringSplitOptions.RemoveEmptyEntries);
        if (splitted.Length >= 2)
        {
            string toSplit = splitted[splitted.Length - 1];
            splitted = toSplit.Split(new[] { '.' }, StringSplitOptions.RemoveEmptyEntries);
            return int.Parse(splitted[0]);
        }

        throw new Exception("Error: nbTrucks could not be read from the file name. Enter it from the command line");
    }
}
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
    private LocalSolver localsolver;

    // Number of customers (= number of nodes minus 1)
    int nbCustomers;

    // Capacity of the trucks
    private int truckCapacity;

    // Demand on each node
    private long[] demands;

    // Distance matrix
    private long[][] distanceMatrix;
    
    // Distances between customers and warehouse
    private long[] distanceWarehouses;

    // Number of trucks
    private int nbTrucks;

    // Decision variables
    private LSExpression[] customersSequences;

    // Are the trucks actually used
    private LSExpression[] trucksUsed;
    
    // Distance travelled by each truck
    LSExpression[] routeDistances;

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

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

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

    // Reads instance data.
    private void readInstance(String fileName) throws IOException {
        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);
    }

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

        trucksUsed = new LSExpression[nbTrucks];
        customersSequences = new LSExpression[nbTrucks];
        routeDistances = new LSExpression[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 the trucks
        model.constraint(model.partition(customersSequences));

        // Create demands and distances as arrays to be able to access it with an "at" operator
        LSExpression demandsArray = model.array(demands);
        LSExpression distanceWarehouseArray = model.array(distanceWarehouses);
        LSExpression distanceArray = model.array(distanceMatrix);

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

            // A truck is used if it visits at least one customer
            trucksUsed[k] = model.gt(c, 0);

            // The quantity needed in each route must not exceed the truck capacity
            LSExpression demandSelector = model.function(i -> model.at(demandsArray, model.at(sequence, i)));
            LSExpression routeQuantity = model.sum(model.range(0, c), demandSelector);
            model.constraint(model.leq(routeQuantity, truckCapacity));

            // Distance traveled by truck k
            LSExpression distSelector = model.function(i -> model.at(
                        distanceArray,
                        model.at(sequence, model.sub(i, 1)),
                        model.at(sequence, i)));
            routeDistances[k] = model.sum(model.sum(model.range(1, c), distSelector),
                                    model.iif(model.gt(c, 0), model.sum(
                                        model.at(distanceWarehouseArray, model.at(sequence, 0)),
                                        model.at(distanceWarehouseArray, model.at(sequence, model.sub(c, 1)))),  0));
        }

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

        // Objective: minimize the number of trucks used, then minimize the distance traveled
        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)
    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-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.print((customersCollection.get(i) + 2) + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "Augerat" format.
    private void readInputCvrp(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("CAPACITY")) {
                    truckCapacity = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("EDGE_WEIGHT_TYPE")) {
                    if (splitted[1].trim().compareTo("EUC_2D") != 0) {
                        throw new RuntimeException("Edge Weight Type " + splitted[1] + " is not supported (only EUC_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 orginal customer indices are in 2..nbNodes 
                    customersX[n - 2] = input.nextInt(); 
                    customersY[n - 2] = input.nextInt(); 
                }    
            }

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

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

            demands = 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("Warehouse demand is supposed to be 0");
                } else {
                    // -2 because orginal customer indices are in 2..nbNodes 
                    demands[n - 2] = demand;
                }
            }

            splitted = 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 warehouseId = input.nextInt();
            if (warehouseId != 1) throw new IOException("Warehouse id is supposed to be 1");

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

    // Computes the distance matrix
    private void computeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY) {
        distanceMatrix = new long[nbCustomers][nbCustomers];
        for (int i = 0; i < nbCustomers; i++) {
            distanceMatrix[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; j++) {
                long dist = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distanceMatrix[i][j] = dist;
                distanceMatrix[j][i] = dist;
            }
        }
        
        distanceWarehouses = new long[nbCustomers];
        for(int i = 0; i < nbCustomers; ++i) {
            distanceWarehouses[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);
    }

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

        throw new RuntimeException("Error: nbTrucks could not be read from the file name. Enter it from the command line");
    }

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

        try {
            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);
            }
        } catch(Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}