Multi Depot Vehicle Routing Problem (MDVRP)

Principles learned

  • Use lists and decision variables

  • Use a lambda expression to define a recursive array

  • Use the partition operator

  • Use the sum operator

Problem

../_images/multi-depot-vehicle-routing-problem-mdvrp.svg

In the Multi Depot Vehicle Routing Problem (MDVRP), which is an extension of the Capacitated Vehicle Routing Problem (CVRP), a fleet of delivery vehicles must serve customers with known demand and service time for a single commodity. Contrary to the CVRP, which treats the MDVRP problem with only one depot location, different depot locations are available, each of them having its own kind of trucks with different properties (demand capacity and service time capacity). The vehicles must end their routes at their starting depot. Each customer can only be served by one vehicle. The objective is to determine which customers should be served by each truck from each known depot by minimizing the total distance traveled by all the trucks, such that all customers are served (demand and service time), the demand for every truck must not exceed its capacity and the time it spends outside its depot must not exceed its time capacity.

Download the example


Data

The instances used in this example come from the Cordeau_2011 instances.

The format of the data files is as follows:

  • One line for: data type (2 for mdvrp), number of trucks per depot, number of customers, number of depots.

  • One line for each depot with: route duration capacity of trucks, quantity capacity of trucks.

  • One line for each customer with: customer id, coordinate x, coordinate y, service time needed, demand needed (we ignore the rest of the line which is for PVRP).

  • One line for each depot with: depot id, coodinate x, coordinate y.

Program

This Hexaly model defines a two dimensional vector customersSequences with size nbTrucks x nbDepots where customersSequences[d][k] is a list which corresponds to the sequence of customers visited by this truck of this depot. To ensure that all customers must be served by only one truck, all the list variables are constrained to form a partition by using the “partition” operator.

The quantity delivered to a customer has to meet the corresponding demand, the total quantity in a truck should not exceed its capacity and the time that the truck spends outside should not exceed its service time capacity.

The quantity served routeQuantity by a truck is defined using a variadic sum operator and the access to the demands array with the “at” operator. We can then constrain this expression to never exceed the capacity of the truck.

In the same way, the two dimensional vector routeDistances and especially routeDistances[d][k] contains the distance traveled by the truck k of depot d. We define it again with a sum accessing the distanceMatrix array with the “at” operator. distanceMatrix gives the distance between two customers, to compute the total distance traveled by a truck we need to add the distances to leave and come back to the depot stored in the distanceDepot array.

If we add to the distance the service times needed to serve every customer this truck should serve, we obtain the time it needs to spend outside (again with the sum and the routeDurationCapacity vector). This one must not exceed its service time capacity.

The objective is to minimize the sum of the elements in routeDistances: totalDistance, because it is the total distance traveled by all the trucks.

Execution:
localsolver mdvrp.lsp inFileName=instances/p01 [lsTimeLimit=] [outputFileName=]
use io;

/* Read instance data */
function input() {
    local usage = "Usage: localsolver mdvrp.lsp inFileName=inputFile" 
            + "[solFileName=outputFile] [lsTimeLimit=timeLimit]";
    if (inFileName == nil) throw usage;

    // Read the data and compute matrices
    readInputMdvrp();
    computeDistanceMatrix();
}

/* Declare the optimization model */
function model() {
    // Sequence of customers visited by each truck
    customersSequences[d in 0...nbDepots][k in 0...nbTrucksPerDepot] <- list(nbCustomers);

    // All customers must be visited by exactly one truck
    constraint partition[d in 0...nbDepots][k in 0...nbTrucksPerDepot](customersSequences[d][k]);

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

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

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

        // We add service Time
        routeServiceTime <- sum(0...c, i => serviceTime[sequence[i]]);
        routeDuration <- routeDistances[d][k] + routeServiceTime;

        // The total distance of this truck should not exceed the duration capacity of the truck
        // (only if we define such a capacity)
        if (routeDurationCapacity[d] > 0) {
            constraint routeDuration <= routeDurationCapacity[d];
        }
    }

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

    // Objective: minimize the total distance traveled
    minimize totalDistance;
}

function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 20;
}

// Write the solution in a file with the following format:
//  - instance, time_limit, total distance
//  - for each depot and for each truck in this depot, the customers visited 
function output() {
    if (solFileName == nil) return;
    local resfile = io.openWrite(solFileName);

    resfile.println("Instance: "+ inFileName + " ; time_limit: " +
            lsTimeLimit + " ; Objective value: " + totalDistance.value);
    for [d in 0...nbDepots] {
        nbTrucksUsed = 0;
        for [k in 0...nbTrucksPerDepot] {
            if (customersSequences[d][k].value.count() > 0) {
                nbTrucksUsed += 1;
            }
        }
        if (nbTrucksUsed > 0){
            n = 1;
            resfile.println("Depot " + (d + 1));
            for [k in 0...nbTrucksPerDepot] {
                if (customersSequences[d][k].value.count() > 0) {
                    resfile.print("Truck " + n + " : ");
                    for [customer in customersSequences[d][k].value] {
                        resfile.print(customer + 1, " ");
                    }
                    n += 1;
                    resfile.println();
                }
            }
            resfile.println();
        }
    }
}

// Input files following "Cordeau"'s format
function readInputMdvrp() {
    local inFile = io.openRead(inFileName);
    inFile.readInt();  // We ignore the first int of the instance

    // Numbers of trucks per depot, customers and depots
    nbTrucksPerDepot = inFile.readInt();
    nbCustomers = inFile.readInt();
    nbDepots = inFile.readInt();

    for [d in 0...nbDepots] {
        routeDurationCapacity[d] = inFile.readInt(); // Time capacity for every type of truck from every depot
        truckCapacity[d] = inFile.readInt(); // Capacity for every type of truck from every depot
    }

    for [n in 0...nbCustomers] {
        inLine = inFile.readln();
        line = inLine.split();

        // Coordinates X and Y, service time and demand for customers
        nodesX[n] = line[1].toDouble();
        nodesY[n] = line[2].toDouble();
        serviceTime[n] = line[3].toInt();
        demands[n] = line[4].toInt();
    }

    for [d in 0...nbDepots] {
        inLine = inFile.readln();
        line = inLine.split();

        // Coordinates X and Y for depots
        depotX[d] = line[1].toDouble();
        depotY[d] = line[2].toDouble();
    }
}

// Compute the distance matrix for customers and for depots/warehouses
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;
        }
    }

    for [i in 0...nbCustomers][d in 0...nbDepots] {
        local localDistance = computeDistDepot(d, i);
        distanceDepot[d][i] = localDistance;
    }
}

// Compute the distance between two points
function computeDist(i, j) {
    return sqrt(pow(nodesX[i] - nodesX[j], 2) + pow(nodesY[i] - nodesY[j], 2));
}

function computeDistDepot(d, j) {
    return sqrt(pow(depotX[d] - nodesX[j], 2) + pow(depotY[d] - nodesY[j], 2));
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python mdvrp.py instances/p01
Execution (Linux)
export PYTHONPATH=/opt/localsolver_13_0/bin/python
python mdvrp.py instances/p01
import localsolver
import sys
import math


def main(instance_file, str_time_limit, output_file):
    #
    # Read instance data
    #
    nb_trucks_per_depot, nb_customers, nb_depots, route_duration_capacity_data, \
        truck_capacity_data, demands_data, service_time_data, \
        distance_matrix_customers_data, distance_warehouse_data = read_input_mdvrp(instance_file)

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

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

        # Vectorization for partition constraint
        customer_sequences_constraint = [customer_sequences[d][k]
                                         for d in range(nb_depots) for k in range(nb_trucks_per_depot)]

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

        # Create LocalSolver arrays to be able to access them with an "at" operator
        demands = model.array(demands_data)
        service_time = model.array(service_time_data)
        dist_customers = model.array(distance_matrix_customers_data)

        # Distances traveled by each truck from each depot
        route_distances = [[None for _ in range(nb_trucks_per_depot)] for _ in range(nb_depots)]

        # Total distance traveled
        total_distance = model.sum()

        for d in range(nb_depots):
            dist_depot = model.array(distance_warehouse_data[d])
            for k in range(nb_trucks_per_depot):
                sequence = customer_sequences[d][k]
                c = model.count(sequence)

                # The quantity needed in each route must not exceed the truck capacity
                demand_lambda = model.lambda_function(lambda j: demands[j])
                route_quantity = model.sum(sequence, demand_lambda)
                model.constraint(route_quantity <= truck_capacity_data[d])

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

                # We add service Time
                service_lambda = model.lambda_function(lambda j: service_time[j])
                route_service_time = model.sum(sequence, service_lambda)

                total_distance.add_operand(route_distances[d][k])

                # The total distance should not exceed the duration capacity of the truck
                # (only if we define such a capacity)
                if (route_duration_capacity_data[d] > 0):
                    model.constraint(route_distances[d][k] + route_service_time <= route_duration_capacity_data[d])

        # Objective: minimize the total distance traveled
        model.minimize(total_distance)

        model.close()

        # Parameterize the solver
        ls.param.time_limit = int(str_time_limit)

        ls.solve()

        #
        # Write the solution in a file with the following format:
        #  - instance, time_limit, total distance
        #  - for each depot and for each truck in this depot, the customers visited
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("Instance: " + instance_file + " ; " + "time_limit: " + str_time_limit + " ; " +
                        "Objective value: " + str(total_distance.value))
                f.write("\n")
                for d in range(nb_depots):
                    trucks_used = []
                    for k in range(nb_trucks_per_depot):
                        if (len(customer_sequences[d][k].value) > 0):
                            trucks_used.append(k)
                    if len(trucks_used) > 0:
                        f.write("Depot " + str(d + 1) + "\n")
                        for k in range(len(trucks_used)):
                            f.write("Truck " + str(k + 1) + " : ")
                            customers_collection = customer_sequences[d][trucks_used[k]].value
                            for p in range(len(customers_collection)):
                                f.write(str(customers_collection[p] + 1) + " ")
                            f.write("\n")
                        f.write("\n")


# Input files following "Cordeau"'s format
def read_input_mdvrp(filename):
    with open(filename) as f:
        instance = f.readlines()

    nb_line = 0
    datas = instance[nb_line].split()

    # Numbers of trucks per depot, customers and depots
    nb_trucks_per_depot = int(datas[1])
    nb_customers = int(datas[2])
    nb_depots = int(datas[3])

    route_duration_capacity = [None]*nb_depots  # Time capacity for every type of truck from every depot
    truck_capacity = [None]*nb_depots  # Capacity for every type of truck from every depot

    for d in range(nb_depots):
        nb_line += 1
        capacities = instance[nb_line].split()

        route_duration_capacity[d] = int(capacities[0])
        truck_capacity[d] = int(capacities[1])

    # Coordinates X and Y, service time and demand for customers
    nodes_xy = [[None, None]] * nb_customers
    service_time = [None] * nb_customers
    demands = [None] * nb_customers

    for n in range(nb_customers):
        nb_line += 1
        customer = instance[nb_line].split()

        nodes_xy[n] = [float(customer[1]), float(customer[2])]

        service_time[n] = int(customer[3])
        demands[n] = int(customer[4])

    # Coordinates X and Y of every depot
    depot_xy = [None] * nb_depots

    for d in range(nb_depots):
        nb_line += 1
        depot = instance[nb_line].split()

        depot_xy[d] = [float(depot[1]), float(depot[2])]

    # Compute the distance matrices
    distance_matrix_customers = compute_distance_matrix_customers(nodes_xy)
    distance_warehouse = compute_distance_warehouse(depot_xy, nodes_xy)

    return nb_trucks_per_depot, nb_customers, nb_depots, route_duration_capacity, \
        truck_capacity, demands, service_time, distance_matrix_customers, distance_warehouse


# Compute the distance matrix for customers
def compute_distance_matrix_customers(nodes_xy):
    nb_customers = len(nodes_xy)
    distance_matrix = [[0 for _ in range(nb_customers)] for _ in range(nb_customers)]
    for i in range(nb_customers):
        for j in range(i+1, nb_customers):
            distij = compute_dist(nodes_xy[i], nodes_xy[j])
            distance_matrix[i][j] = distij
            distance_matrix[j][i] = distij
    return distance_matrix


# Compute the distance matrix for warehouses/depots
def compute_distance_warehouse(depot_xy, nodes_xy):
    nb_customers = len(nodes_xy)
    nb_depots = len(depot_xy)
    distance_warehouse = [[0 for _ in range(nb_customers)] for _ in range(nb_depots)]

    for i in range(nb_customers):
        for d in range(nb_depots):
            distance_warehouse[d][i] = compute_dist(depot_xy[d], nodes_xy[i])

    return distance_warehouse


# Compute the distance between two points
def compute_dist(p, q):
    return math.sqrt(math.pow(p[0] - q[0], 2) + math.pow(p[1] - q[1], 2))


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python mdvrp.py input_file [output_file] [time_limit]")
        sys.exit(1)

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

    main(instance_file, str_time_limit, output_file)
Compilation / Execution (Windows)
cl /EHsc mdvrp.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver130.lib
mdvrp instances/p01.dat
Compilation / Execution (Linux)
g++ mdvrp.cpp -I/opt/localsolver_13_0/include -llocalsolver130 -lpthread -o mdvrp
./mdvrp instances/p01
#include "localsolver.h"

#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace localsolver;
using namespace std;

class Mdvrp {
public:
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    // Number of trucks per depot
    int nbTrucksPerDepot;

    // Capacity of the trucks per depot
    vector<int> truckCapacity;

    // Duration capacity of the trucks per depot
    vector<int> routeDurationCapacity;

    // Service time per customer
    vector<int> serviceTimeData;

    // Demand per customer
    vector<int> demandsData;

    // Distance matrix between customers
    vector<vector<double>> distanceMatrixCustomersData;

    // Distances between customers and depots
    vector<vector<double>> distanceWarehouseData;

    // Decision variables
    vector<vector<LSExpression>> customersSequences;

    // Distance traveled by all the trucks
    LSExpression totalDistance;

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

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

        // Sequence of customers visited by each truck
        customersSequences.resize(nbDepots);

        // Vectorization for partition constraint
        vector<LSExpression> customersSequencesConstraint(nbDepots * nbTrucksPerDepot);

        for (int d = 0; d < nbDepots; ++d) {
            customersSequences[d].resize(nbTrucksPerDepot);
            for (int k = 0; k < nbTrucksPerDepot; ++k) {
                customersSequences[d][k] = model.listVar(nbCustomers);
                customersSequencesConstraint[d * nbTrucksPerDepot + k] = customersSequences[d][k];
            }
        }

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

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = model.array(demandsData.begin(), demandsData.end());
        LSExpression serviceTime = model.array(serviceTimeData.begin(), serviceTimeData.end());

        LSExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            distMatrix.addOperand(
                model.array(distanceMatrixCustomersData[n].begin(), distanceMatrixCustomersData[n].end()));
        }

        // Distances traveled by each truck from each depot
        vector<vector<LSExpression>> routeDistances;
        routeDistances.resize(nbDepots);

        // Total distance traveled
        totalDistance = model.sum();

        for (int d = 0; d < nbDepots; ++d) {
            routeDistances[d].resize(nbTrucksPerDepot);
            LSExpression distDepot = model.array(distanceWarehouseData[d].begin(), distanceWarehouseData[d].end());

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

                // The quantity needed in each route must not exceed the truck capacity
                LSExpression demandLambda = model.createLambdaFunction([&](LSExpression j) { return demands[j]; });

                LSExpression routeQuantity = model.sum(sequence, demandLambda);

                model.constraint(routeQuantity <= truckCapacity[d]);

                // Distance traveled by truck k of depot d
                LSExpression distLambda = model.createLambdaFunction(
                    [&](LSExpression i) { return model.at(distMatrix, sequence[i - 1], sequence[i]); });

                routeDistances[d][k] = model.sum(model.range(1, c), distLambda) +
                                       model.iif(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

                totalDistance.addOperand(routeDistances[d][k]);

                // We add service time
                LSExpression serviceLambda = model.createLambdaFunction([&](LSExpression j) { return serviceTime[j]; });
                LSExpression routeServiceTime = model.sum(sequence, serviceLambda);

                // The total distance should not exceed the duration capacity of the truck
                // (only if we define such a capacity)
                if (routeDurationCapacity[d] > 0) {
                    model.constraint(routeDistances[d][k] + routeServiceTime <= routeDurationCapacity[d]);
                }
            }
        }

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

        model.close();

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    // Write the solution in a file with the following format:
    //  - instance, time_limit, total distance
    //  - for each depot and for each truck in this depot, the customers visited
    void writeSolution(const string& fileName, const string& instanceFile, const string& timeLimit) {
        ofstream outfile;
        outfile.open(fileName.c_str());

        outfile << "Instance: " << instanceFile << " ; time_limit: " << timeLimit
                << " ; Objective value: " << totalDistance.getDoubleValue() << endl;
        for (int d = 0; d < nbDepots; ++d) {
            vector<int> trucksUsed;
            for (int k = 0; k < nbTrucksPerDepot; ++k) {
                if (customersSequences[d][k].getCollectionValue().count() > 0) {
                    trucksUsed.push_back(k);
                }
            }
            if (trucksUsed.size() > 0) {
                outfile << "Depot " << d + 1 << endl;
                for (int k = 0; k < trucksUsed.size(); ++k) {
                    outfile << "Truck " << (k + 1) << " : ";
                    LSCollection customersCollection = customersSequences[d][trucksUsed[k]].getCollectionValue();
                    for (int p = 0; p < customersCollection.count(); ++p) {
                        outfile << customersCollection[p] + 1 << " ";
                    }
                    outfile << endl;
                }
                outfile << endl;
            }
        }
    }

private:
    // Input files following "Cordeau"'s format
    void readInputMdvrp(const string& fileName) {
        ifstream infile;
        infile.open(fileName.c_str());

        infile.ignore(); // We ignore the first int of the instance

        // Numbers of trucks per depot, customers and depots
        infile >> nbTrucksPerDepot;
        infile >> nbCustomers;
        infile >> nbDepots;

        routeDurationCapacity.resize(nbDepots);
        truckCapacity.resize(nbDepots);

        for (int d = 0; d < nbDepots; ++d) {
            infile >> routeDurationCapacity[d];
            infile >> truckCapacity[d];
        }

        // Coordinates X and Y, service time and demand for customers
        vector<double> nodesX(nbCustomers);
        vector<double> nodesY(nbCustomers);
        serviceTimeData.resize(nbCustomers);
        demandsData.resize(nbCustomers);

        for (int n = 0; n < nbCustomers; ++n) {
            int id;
            int bin;
            infile >> id;
            infile >> nodesX[id - 1];
            infile >> nodesY[id - 1];
            infile >> serviceTimeData[id - 1];
            infile >> demandsData[id - 1];

            // Ignore the end of the line
            infile.ignore(numeric_limits<streamsize>::max(), '\n');
        }
        // Coordinates X and Y for depots
        vector<double> depotX(nbDepots);
        vector<double> depotY(nbDepots);

        for (int d = 0; d < nbDepots; ++d) {
            int id;
            int bin;
            infile >> id;
            infile >> depotX[id - nbCustomers - 1];
            infile >> depotY[id - nbCustomers - 1];

            // Ignore the end of the line
            infile.ignore(numeric_limits<streamsize>::max(), '\n');
        }

        // Compute the distance matrices
        computeDistanceMatrixCustomers(nodesX, nodesY);
        computeDistanceWarehouse(depotX, depotY, nodesX, nodesY);

        infile.close();
    }

    // Compute the distance matrix for customers
    void computeDistanceMatrixCustomers(const vector<double>& nodesX, const vector<double>& nodesY) {
        distanceMatrixCustomersData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distanceMatrixCustomersData[i].resize(nbCustomers);
        }

        for (int i = 0; i < nbCustomers; ++i) {
            distanceMatrixCustomersData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                double distance = computeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrixCustomersData[i][j] = distance;
                distanceMatrixCustomersData[j][i] = distance;
            }
        }
    }

    // Compute the distance matrix for warehouses
    void computeDistanceWarehouse(const vector<double>& depotX, const vector<double>& depotY,
                                  const vector<double>& nodesX, const vector<double>& nodesY) {
        distanceWarehouseData.resize(nbDepots);

        for (int d = 0; d < nbDepots; ++d) {
            distanceWarehouseData[d].resize(nbCustomers);
            for (int i = 0; i < nbCustomers; ++i) {
                distanceWarehouseData[d][i] = computeDist(nodesX[i], depotX[d], nodesY[i], depotY[d]);
            }
        }
    }

    // Compute the distance between two points
    double computeDist(double xi, double xj, double yi, double yj) { return sqrt(pow(xi - xj, 2) + pow(yi - yj, 2)); }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: mdvrp inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }
    const char* instanceFile = argv[1];
    const char* outputFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "20";

    try {
        Mdvrp model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));

        // If we want to write the solution
        if (outputFile != NULL)
            model.writeSolution(outputFile, instanceFile, strTimeLimit);
    } catch (const exception& e) {
        cerr << "An error occured: " << e.what() << endl;
    }

    return 0;
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc Mdvrp.cs /reference:localsolvernet.dll
mdvrp instances\p01
using localsolver;
using System;
using System.IO;
using System.Collections.Generic;

public class Mdvrp : IDisposable
{
    // LocalSolver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    //Number of trucks per depot
    int nbTrucksPerDepot;

    // Capacity of the trucks per depot
    long[] truckCapacity;

    // Duration capacity of the trucks per depot
    long[] routeDurationCapacity;

    // Service time per customer
    long[] serviceTimeData;

    // Demand per customer
    long[] demandsData;

    // Distance matrix between customers
    double[][] distanceMatrixCustomersData;

    // Distances between customers and depots
    double[][] distanceWarehouseData;

    // Decision variable
    LSExpression[][] customersSequences;

    // Distance traveled by all the trucks
    LSExpression totalDistance;

    public Mdvrp()
    {
        localsolver = new LocalSolver();
    }

    void ReadInstance(string fileName)
    {
        readInputMdvrp(fileName);
    }

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

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

        // Sequence of customers visited by each truck
        customersSequences = new LSExpression[nbDepots][];

        // Vectorization for partition constraint
        LSExpression[] customersSequencesConstraint = new LSExpression[nbTrucksPerDepot * nbDepots];

        for (int d = 0; d < nbDepots; ++d)
        {
            customersSequences[d] = new LSExpression[nbTrucksPerDepot];
            for (int k = 0; k < nbTrucksPerDepot; ++k)
            {
                customersSequences[d][k] = model.List(nbCustomers);
                customersSequencesConstraint[d * nbTrucksPerDepot + k] = customersSequences[d][k];
            }
        }

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

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = model.Array(demandsData);
        LSExpression serviceTime = model.Array(serviceTimeData);

        LSExpression distMatrix = model.Array(distanceMatrixCustomersData);

        //Distances for each truck from each depot 
        LSExpression[][] routeDistances = new LSExpression[nbDepots][];

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

        for (int d = 0; d < nbDepots; ++d)
        {
            routeDistances[d] = new LSExpression[nbTrucksPerDepot];
            LSExpression distDepot = model.Array(distanceWarehouseData[d]);
            for (int k = 0; k < nbTrucksPerDepot; ++k)
            {
                LSExpression sequence = customersSequences[d][k];
                LSExpression c = model.Count(sequence);

                // The quantity needed in each route must not exceed the truck capacity
                LSExpression demandLambda = model.LambdaFunction(j => demands[j]);
                LSExpression routeQuantity = model.Sum(sequence, demandLambda);
                model.Constraint(routeQuantity <= truckCapacity[d]);

                // Distance traveled by truck k of depot d
                LSExpression distLambda = model.LambdaFunction(
                    i => distMatrix[sequence[i - 1], sequence[i]]
                );
                routeDistances[d][k] =
                    model.Sum(model.Range(1, c), distLambda)
                    + model.If(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

                totalDistance.AddOperand(routeDistances[d][k]);

                // We add service time
                LSExpression serviceLambda = model.LambdaFunction(j => serviceTime[j]);
                LSExpression routeServiceTime = model.Sum(sequence, serviceLambda);

                // The total distance should not exceed the duration capacity of the truck
                // (only if we define such a capacity)
                if (routeDurationCapacity[d] > 0)
                {
                    model.Constraint(routeDistances[d][k] + routeServiceTime <= routeDurationCapacity[d]);
                }
            }
        }

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

        model.Close();

        // Parametrize the solver
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.Solve();
    }

    // Write the solution in a file with the following format:
    //  - instance, time_limit, total distance
    //  - for each depot and for each truck in this depot, the customers visited 
    void WriteSolution(string fileName, string instanceFile, string timeLimit)
    {
        using (StreamWriter outfile = new StreamWriter(fileName))
        {
            outfile.WriteLine("Instance: " + instanceFile + " ; time_limit: " + timeLimit +
                    " ; Objective value: " + totalDistance.GetDoubleValue());
            for (int d = 0; d < nbDepots; ++d)
            {
                List<int> trucksUsed = new List<int>();
                for (int k = 0; k < nbTrucksPerDepot; ++k)
                {
                    if (customersSequences[d][k].GetCollectionValue().Count() > 0)
                    {
                        trucksUsed.Add(k);
                    }
                }

                if (trucksUsed.Count > 0)
                {
                    outfile.WriteLine("Depot " + (d + 1));
                    for (int k = 0; k < trucksUsed.Count; ++k)
                    {
                        outfile.Write("Truck " + (k + 1) + " : ");
                        LSCollection customersCollection = customersSequences[d][trucksUsed[k]].GetCollectionValue();
                        for (int p = 0; p < customersCollection.Count(); ++p)
                            outfile.Write((customersCollection[p] + 1) + " ");
                        outfile.WriteLine();
                    }
                    outfile.WriteLine();
                }
            }
        }
    }

    // Input files following "Cordeau"'s format
    private void readInputMdvrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);

            // Numbers of trucks per depot, customers and depots
            nbTrucksPerDepot = int.Parse(splitted[1]);
            nbCustomers = int.Parse(splitted[2]);
            nbDepots = int.Parse(splitted[3]);

            routeDurationCapacity = new long[nbDepots];
            truckCapacity = new long[nbDepots];

            for (int d = 0; d < nbDepots; ++d)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);

                routeDurationCapacity[d] = int.Parse(splitted[0]);
                truckCapacity[d] = int.Parse(splitted[1]);
            }

            // Coordinates X and Y, service time and demand for customers
            double[] nodesX = new double[nbCustomers];
            double[] nodesY = new double[nbCustomers];
            serviceTimeData = new long[nbCustomers];
            demandsData = new long[nbCustomers];

            for (int n = 0; n < nbCustomers; ++n)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);

                nodesX[n] = double.Parse(splitted[1]);
                nodesY[n] = double.Parse(splitted[2]);
                serviceTimeData[n] = int.Parse(splitted[3]);
                demandsData[n] = int.Parse(splitted[4]);
            }

            // Coordinates X and Y for depots
            double[] DepotX = new double[nbDepots];
            double[] DepotY = new double[nbDepots];

            for (int d = 0; d < nbDepots; ++d)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);

                DepotX[d] = double.Parse(splitted[1]);
                DepotY[d] = double.Parse(splitted[2]);
            }

            // Compute the distance matrices
            ComputeDistanceMatrixCustomers(nodesX, nodesY);
            ComputeDistanceWarehouse(DepotX, DepotY, nodesX, nodesY);
        }
    }

    // Compute the distance matrix for customers
    private void ComputeDistanceMatrixCustomers(double[] nodesX, double[] nodesY)
    {
        distanceMatrixCustomersData = new double[nbCustomers][];

        for (int i = 0; i < nbCustomers; ++i)
            distanceMatrixCustomersData[i] = new double[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distanceMatrixCustomersData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                double dist = ComputeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrixCustomersData[i][j] = dist;
                distanceMatrixCustomersData[j][i] = dist;
            }
        }

    }

    // Compute the distance matrix for depots/warehouses
    private void ComputeDistanceWarehouse(double[] DepotX, double[] DepotY, double[] nodesX, double[] nodesY)
    {
        distanceWarehouseData = new double[nbDepots][];

        for (int d = 0; d < nbDepots; ++d)
        {
            distanceWarehouseData[d] = new double[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
            {
                distanceWarehouseData[d][i] = ComputeDist(nodesX[i], DepotX[d], nodesY[i], DepotY[d]);
            }
        }
    }

    private double ComputeDist(double xi, double xj, double yi, double yj)
    {
        return Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Mdvrp inputFile [outputFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "20";

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

}
Compilation / Execution (Windows)
javac Mdvrp.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. mdvrp instances\p01
Compilation / Execution (Linux)
javac Mdvrp.java -cp /opt/localsolver_13_0/bin/localsolver.jar
java -cp /opt/localsolver_13_0/bin/localsolver.jar:. mdvrp instances/p01
import java.util.*;
import java.io.*;
import localsolver.*;

public class Mdvrp {
    // LocalSolver
    private final LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    // Number of trucks per depot
    private int nbTrucksPerDepot;

    // Capacity of the trucks per depot
    private long[] truckCapacity;

    // Duration capacity of the trucks per depot
    private long[] routeDurationCapacity;

    // Service time per customer
    private long[] serviceTimeData;

    // Demand per customer
    private long[] demandsData;

    // Distance matrix between customers
    private double[][] distanceMatrixCustomersData;

    // Distance matrix between customers and depots
    private double[][] distanceWarehouseData;

    // Decision variable
    private LSExpression[][] customersSequences;

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

    private Mdvrp(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    private void readInstance(String fileName) throws IOException {
        readInputMdvrp(fileName);
    }

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

        // Sequence of customers visited by each truck
        customersSequences = new LSExpression[nbDepots][nbTrucksPerDepot];

        // Vectorization for partition constraint
        LSExpression[] customersSequencesConstraint = new LSExpression[nbDepots * nbTrucksPerDepot];

        // Sequence of customers visited by each truck
        for (int d = 0; d < nbDepots; ++d) {
            for (int k = 0; k < nbTrucksPerDepot; ++k) {
                customersSequences[d][k] = model.listVar(nbCustomers);
                customersSequencesConstraint[d * nbTrucksPerDepot + k] = customersSequences[d][k];
            }
        }
        // All customers must be visited by exactly one truck
        model.constraint(model.partition(customersSequencesConstraint));

        // Create LocalSolver arrays to be able to access them with an "at" operator
        LSExpression demands = model.array(demandsData);
        LSExpression serviceTime = model.array(serviceTimeData);

        LSExpression distMatrix = model.array(distanceMatrixCustomersData);

        // Distances for each truck from each depot
        LSExpression[][] routeDistances = new LSExpression[nbDepots][nbTrucksPerDepot];

        // Total distance traveled
        totalDistance = model.sum();

        for (int d = 0; d < nbDepots; ++d) {
            LSExpression distDepot = model.array(distanceWarehouseData[d]);

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

                // The quantity needed in each route must not exceed the truck capacity
                LSExpression demandLambda = model.lambdaFunction(j -> model.at(demands, j));
                LSExpression routeQuantity = model.sum(sequence, demandLambda);
                model.constraint(model.leq(routeQuantity, truckCapacity[d]));

                // Distance traveled by truck k of depot d
                LSExpression distLambda = model
                        .lambdaFunction(
                                i -> model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)));
                routeDistances[d][k] = model.sum(model.sum(model.range(1, c), distLambda),
                        model.iif(model.gt(c, 0), model.sum(model.at(distDepot, model.at(sequence, 0)),
                                model.at(distDepot, model.at(sequence, model.sub(c, 1)))), 0));

                totalDistance.addOperand(routeDistances[d][k]);

                // We add service time
                LSExpression serviceLambda = model.lambdaFunction(j -> model.at(serviceTime, j));
                LSExpression routeServiceTime = model.sum(sequence, serviceLambda);

                // The total distance should not exceed the duration capacity of the truck (only
                // if we define such a capacity)
                if (routeDurationCapacity[d] > 0) {
                    model.constraint(
                            model.leq(model.sum(routeDistances[d][k], routeServiceTime), routeDurationCapacity[d]));
                }
            }
        }

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

        model.close();

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    // Write the solution in a file with the following format:
    // - instance, time_limit, total distance
    // - for each depot and for each truck in this depot, the customers visited
    private void writeSolution(String fileName, String instanceFile, String timeLimit) throws IOException {
        try (PrintWriter outfile = new PrintWriter(fileName)) {
            outfile.println("Instance: " + instanceFile + " ; time_limit: " + timeLimit + " ; Objective value: "
                    + totalDistance.getDoubleValue());
            for (int d = 0; d < nbDepots; ++d) {
                ArrayList<Integer> trucksUsed = new ArrayList<Integer>();
                for (int k = 0; k < nbTrucksPerDepot; k++) {
                    if (customersSequences[d][k].getCollectionValue().count() > 0) {
                        trucksUsed.add(k);
                    }
                }

                if (trucksUsed.size() > 0) {
                    outfile.println("Depot " + (d + 1));
                    for (int k = 0; k < trucksUsed.size(); ++k) {
                        outfile.print("Truck " + (k + 1) + " : ");
                        LSCollection customersCollection = customersSequences[d][trucksUsed.get(k)]
                                .getCollectionValue();
                        for (int p = 0; p < customersCollection.count(); ++p)
                            outfile.print((customersCollection.get(p) + 1) + " ");
                        outfile.println();
                    }
                    outfile.println();
                }
            }
        }
    }

    // Input files following "Cordeau"'s format
    private void readInputMdvrp(String fileName) throws IOException {

        try (Scanner input = new Scanner(new File(fileName))) {
            String[] splitted;

            splitted = input.nextLine().split(" ");

            // Numbers of trucks per depot, customers and depots
            nbTrucksPerDepot = Integer.parseInt(splitted[1].trim());
            nbCustomers = Integer.parseInt(splitted[2].trim());
            nbDepots = Integer.parseInt(splitted[3].trim());

            routeDurationCapacity = new long[nbDepots];
            truckCapacity = new long[nbDepots];

            for (int d = 0; d < nbDepots; ++d) {
                routeDurationCapacity[d] = input.nextInt();
                truckCapacity[d] = input.nextInt();
            }

            // Coordinates X and Y, service time and demand for customers
            double[] nodesX = new double[nbCustomers];
            double[] nodesY = new double[nbCustomers];
            serviceTimeData = new long[nbCustomers];
            demandsData = new long[nbCustomers];

            for (int n = 0; n < nbCustomers; ++n) {
                input.nextLine();
                input.nextInt();

                nodesX[n] = input.nextDouble();
                nodesY[n] = input.nextDouble();
                serviceTimeData[n] = input.nextInt();
                demandsData[n] = input.nextInt();
            }

            // Coordinates X and Y for depots
            double[] DepotX = new double[nbDepots];
            double[] DepotY = new double[nbDepots];

            for (int d = 0; d < nbDepots; ++d) {
                input.nextLine();
                input.nextInt();
                DepotX[d] = input.nextDouble();
                DepotY[d] = input.nextDouble();
            }

            // Compute the distance matrices
            ComputeDistanceMatrixCustomers(nodesX, nodesY);
            ComputeDistanceWarehouse(DepotX, DepotY, nodesX, nodesY);
        }
    }

    // Compute the distance matrix for customers
    private void ComputeDistanceMatrixCustomers(double[] nodesX, double[] nodesY) {
        distanceMatrixCustomersData = new double[nbCustomers][nbCustomers];

        for (int i = 0; i < nbCustomers; ++i) {
            distanceMatrixCustomersData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                double dist = computeDist(nodesX[i], nodesX[j], nodesY[i], nodesY[j]);
                distanceMatrixCustomersData[i][j] = dist;
                distanceMatrixCustomersData[j][i] = dist;
            }
        }
    }

    // Compute the distance matrix for depots/warehouses
    private void ComputeDistanceWarehouse(double[] DepotX, double[] DepotY, double[] nodesX, double[] nodesY) {
        distanceWarehouseData = new double[nbDepots][nbCustomers];

        for (int d = 0; d < nbDepots; ++d) {
            for (int i = 0; i < nbCustomers; ++i) {
                distanceWarehouseData[d][i] = computeDist(nodesX[i], DepotX[d], nodesY[i], DepotY[d]);
            }
        }
    }

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

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java Mdvrp inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }
        try (LocalSolver localsolver = new LocalSolver()) {
            String instanceFile = args[0];
            String outputFile = args.length > 1 ? args[1] : null;
            String strTimeLimit = args.length > 2 ? args[2] : "20";

            Mdvrp model = new Mdvrp(localsolver);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));

            if (outputFile != null) {
                model.writeSolution(outputFile, instanceFile, strTimeLimit);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }

}