Location Routing Problem (LRP)

Principles learned

  • Combine list and set decision variables

  • Use the find operator

  • Use a lambda expression to define a recursive array

Problem

../_images/location-routing-problem-lrp.svg

In the Location Routing Problem (LRP), which is an extension of the Capacitated Vehicle Routing Problem (CVRP), a fleet of delivery vehicles with uniform capacity must service customers with known demand for a single commodity. Contrary to the CVRP, which treats the LRP problem with only one depot location, different depot locations are available, each of them having its own opening cost and capacity. The vehicles end their routes at their starting depot. Each customer can only be served by one vehicle. The objectives are to determine which depots to open and routes to take while minimizing the total cost, such that all customers are served, and the total demand served by each depot and truck does not exceed its capacity.

Download the example


Data

The instances used in this example come from the S. Barreto instances.

The format of the data files is as follows:

  • The number of customers

  • The number of depots

  • The x and y coordinates of the depots and the customers

  • The capacity of the delivery vehicles

  • The capacity of each depot

  • The demand for every customer

  • The opening cost of each depot

  • The opening cost of a route

  • The areCostDouble boolean value indicating if costs should be rounded

Program

This LocalSolver model defines a sequence for nbTrucks trucks as a list variable (customersSequences[r]). 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 by using to the “partition” operator. The model also defines a group of sequences associated to every depot as a set variable (depots[d]). The depots are also constrained to form a partition, to ensure that every sequence is associated to a depot.

The number of sequences nbTrucks is determined as the minimum number of trucks necessary to serve every customer without overloading the vehicles, multiplied by a factor 1.5 to ensure feasibility. The quantity delivered to a customer has to meet the corresponding demand. The quantity served in a sequence quantityServed[r] is defined with a sum and the access to the demands array with the “at” operator and must not exceed its capacity.

To determine the depot assigned to a route, we use the “find” operator. A depot is considered opened if at least one (non-empty) sequence uses it. This condition is computed with the “count” operator. The depotsCost are defined as the sum of the opening costs of the opened depots. Similarly, a sequence is considered opened if at least one customer is served for this sequence. We calculate the total distance for every sequence by adding the distance from the depot to the first customer, the distances between every customer and the distance back to the depot from the last customer. The sum of sequenceLength and openingCostRoute corresponds to the routingCost.

The objective is to minimize the sum of depotsCost and routingCost.

Execution:
localsolver location_routing_problem.lsp inFileName=instances/coordChrist100.dat [lsTimeLimit=] [solFileName=]
/********** location_routing_problem.lsp **********/

use io;

function input() {

    local usage = "Usage : localsolver location_routing_problem.lsp inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]";
    if (inFileName == nil) throw usage;

    readInputLrp();
    computeDistances();

    minNbTrucks = ceil(sum[c in 0..nbCustomers-1](demands[c]) / truckCapacity);
    nbTrucks = ceil(1.5 * minNbTrucks);
}

function model() {
    // A route is represented as a list containing the customers in the order they are visited
    customersSequences[r in 0..nbTrucks-1] <- list(nbCustomers);   
    // A depot is represented as a set containing the associated customersSequences
    depots[d in 0..nbDepots-1] <- set(nbTrucks);  
    
    // All customers should be assigned to a route
    constraint partition(customersSequences);  
    // All the customersSequences should be assigned to a depot
    constraint partition(depots);  
    
    for [r in 0..nbTrucks-1] {
        local sequence <- customersSequences[r];
        local c <- count(sequence);

        quantityServed[r] <- sum(sequence, c => demands[c]); 
        // The quantity needed in each route must not exceed the vehicle capacity
        constraint quantityServed[r] <= truckCapacity;
        // A route is used if it serves at least one customer  
        sequenceUsed[r] <- c > 0;

        // The "find" function gets the depot that is assigned to the route 
        associatedDepot[r] <- find(depots, r);  
        
        sequenceLength[r] <- sum(0..c - 2, i => distCustomers[sequence[i]][sequence[i + 1]])
                + (sequenceUsed[r] ? distCustomersDepots[sequence[0]][associatedDepot[r]] 
                + distCustomersDepots[sequence[c - 1]][associatedDepot[r]] : 0);
        // The route cost is the sum of the opening cost and the route length
        routeCost[r] <- sequenceLength[r] + sequenceUsed[r] * openingRouteCost;
    }

    for [d in 0..nbDepots-1]{
        // The total demand served by a depot must not exceed its capacity
        constraint sum(depots[d], i => quantityServed[i]) <= depotsCapacity[d];
        // A depot is open if at least a route starts from there
        isDepotOpen[d] <- count(depots[d]) > 0;  
    }

    depotsCost <- sum[d in 0..nbDepots-1](isDepotOpen[d] * openingCostDepots[d]);
    routingCost <- sum[r in 0..nbTrucks-1](routeCost[r]);

    totalCost <- routingCost + depotsCost;
    minimize totalCost;
}

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

function output() {

    if (solFileName == nil) return;
    local resFile = io.openWrite(solFileName);
    resFile.println("File name: ", inFileName + "; totalCost = " + totalCost.value);
    for [r in 0..nbTrucks-1]{
        if (sequenceUsed[r].value){
            resFile.print("Route ", r, ", assigned to depot ", associatedDepot[r].value, " : ");
            for [customer in customersSequences[r].value]{
                resFile.print(customer, " ");
            }
            resFile.println();
        }
    }
    resFile.close();
}

function readInputLrp() {
    if (inFileName.endsWith(".dat")) readInputDat();
    else throw("Unknown file format");
}

function readInputDat() {
    local inFile = io.openRead(inFileName);

    nbCustomers = inFile.readInt();
    nbDepots = inFile.readInt();

    coordinatesDepots[0..nbDepots-1][0..1] = inFile.readDouble();
    coordinatesCustomers[0..nbCustomers-1][0..1] = inFile.readDouble();

    truckCapacity = inFile.readInt();
    depotsCapacity[0..nbDepots-1] = inFile.readInt();

    demands[0..nbCustomers-1] = inFile.readInt();

    local tempOpeningCostDepots[0..nbDepots-1] = inFile.readDouble();
    local tempOpeningCostRoute = inFile.readDouble();

    areCostsDouble = inFile.readInt();

    if (areCostsDouble == 1) {
        openingCostDepots[d in 0..nbDepots-1] = tempOpeningCostDepots[d];
        openingRouteCost = tempOpeningCostRoute;
    } else {
        openingCostDepots[d in 0..nbDepots-1] = round(tempOpeningCostDepots[d]);
        openingRouteCost = round(tempOpeningCostRoute);
    }
}



function computeDistances(){
    local tempDistanceCustomers[c0 in 0..nbCustomers-1][c1 in 0..nbCustomers-1] = 
            sqrt(pow((coordinatesCustomers[c0][0] - coordinatesCustomers[c1][0]), 2)
                    + pow((coordinatesCustomers[c0][1] - coordinatesCustomers[c1][1]), 2));
    local tempDistanceCustomersDepots[c in 0..nbCustomers-1][d in 0..nbDepots-1] = 
            sqrt(pow((coordinatesCustomers[c][0] - coordinatesDepots[d][0]), 2)
                    + pow((coordinatesCustomers[c][1] - coordinatesDepots[d][1]), 2));

    if (areCostsDouble == 1) {
        distCustomers[c0 in 0..nbCustomers-1][c1 in 0..nbCustomers-1] = tempDistanceCustomers[c0][c1];
        distCustomersDepots[c in 0..nbCustomers-1][d in 0..nbDepots-1] = tempDistanceCustomersDepots[c][d];
    } else {
        distCustomers[c0 in 0..nbCustomers-1][c1 in 0..nbCustomers-1] = ceil(100 * tempDistanceCustomers[c0][c1]);
        distCustomersDepots[c in 0..nbCustomers-1][d in 0..nbDepots-1] = ceil(100 * tempDistanceCustomersDepots[c][d]);
    }
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python location_routing_problem.py instances/coordChrist100.dat
Execution (Linux)
export PYTHONPATH=/opt/localsolver_11_0/bin/python
python location_routing_problem.py instances/coordChrist100.dat
########## location_routing_problem.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):
    #
    # Reads instance data
    #
    (nb_customers, nb_depots, vehicle_capacity, opening_route_cost, demands, capacity_depots,
     opening_depots_cost, dist_customers, dist_customers_warehouses) = read_input_lrp(instance_file)

    min_nb_trucks = int(math.ceil(sum(demands) / vehicle_capacity))
    nb_trucks = int(math.ceil(1.5 * min_nb_trucks))

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

        # A sequence is represented as a list containing the customers in the order they are visited
        customers_sequences = [m.list(nb_customers) for _ in range(nb_trucks)]
        # All customers should be affected to a sequence
        m.constraint(m.partition(customers_sequences))

        # A depot is represented as a set containing the associated sequences
        depots = [m.set(nb_trucks) for _ in range(nb_depots)]
        # All the sequences should be affected to a depot
        m.constraint(m.partition(depots))

        quantity_served = [None] * nb_trucks
        sequence_cost = [None] * nb_trucks
        sequence_used = [None] * nb_trucks
        sequence_length = [None] * nb_trucks
        associated_depot = [None] * nb_trucks

        demands = m.array(demands)
        dist_customers_array = m.array()
        dist_customers_depots_array = m.array()

        # We copy the python tables as LS Arrays
        for i in range(nb_customers):
            dist_customers_array.add_operand(m.array(dist_customers[i]))
            dist_customers_depots_array.add_operand(m.array(dist_customers_warehouses[i]))

        for r in range(nb_trucks):
            sequence = customers_sequences[r]
            c = m.count(sequence)

            # A sequence is used if it serves at least one customer
            sequence_used[r] = c > 0
            # The "find" function gets the depot that is assigned to the sequence
            associated_depot[r] = m.find(m.array(depots), r)

            # The quantity needed in each sequence must not exceed the vehicle capacity
            demand_selector = m.lambda_function(lambda i: demands[i])
            quantity_served[r] = m.sum(sequence, demand_selector)
            m.constraint(quantity_served[r] <= vehicle_capacity)

            # Distance traveled by each truck
            dist_selector = m.lambda_function(lambda i: m.at(dist_customers_array, sequence[i], sequence[i + 1]))
            depot = associated_depot[r]
            sequence_length[r] = m.sum(m.range(0, c - 1), dist_selector) + \
                m.iif(sequence_used[r],
                      m.at(dist_customers_depots_array, sequence[0], depot) +
                          m.at(dist_customers_depots_array, sequence[c-1], depot),
                      0)

            # The sequence cost is the sum of the opening cost and the sequence length
            sequence_cost[r] = sequence_used[r] * opening_route_cost + sequence_length[r]

        depot_cost = [None] * nb_depots
        quantity_served_array = m.array(quantity_served)
        for d in range(nb_depots):
            # A depot is open if at least one sequence starts from there
            depot_cost[d] = (m.count(depots[d]) > 0) * opening_depots_cost[d]

            # The total demand served by a depot must not exceed its capacity
            depot_select = m.lambda_function(lambda r: quantity_served_array[r])
            depot_quantity = m.sum(depots[d], depot_select)
            m.constraint(depot_quantity <= capacity_depots[d])

        depots_cost = m.sum(depot_cost)
        routing_cost = m.sum(sequence_cost)
        totalCost = routing_cost + depots_cost

        m.minimize(totalCost)

        m.close()

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

        ls.solve()

        if (sol_file != None):
            with open(sol_file, 'w') as file:
                file.write("File name: %s; totalCost = %d \n" % (instance_file, totalCost.value))
                for r in range(nb_trucks):
                    if (sequence_used[r].value):
                        file.write("Route %d, assigned to depot %d: " % (r, associated_depot[r].value))
                        for customer in customers_sequences[r].value:
                            file.write("%d " % customer)
                        file.write("\n")


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

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

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

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

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

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

    temp_opening_cost_depot = [None] * nb_depots
    for i in range(nb_depots):
        temp_opening_cost_depot[i] = float(next(file_it))
    temp_opening_route_cost = int(next(file_it))
    are_cost_double = int(next(file_it))

    opening_depots_cost = [None] * nb_depots
    if are_cost_double == 1:
        opening_depots_cost = temp_opening_cost_depot
        opening_route_cost = temp_opening_route_cost
    else:
        opening_route_cost = round(temp_opening_route_cost)
        for i in range(nb_depots):
            opening_depots_cost[i] = round(temp_opening_cost_depot[i])

    distance_customers = compute_distance_matrix(x_customer, y_customer, are_cost_double)
    distance_customers_depots = compute_distance_depot(x_customer, y_customer, x_depot, y_depot, are_cost_double)

    return (nb_customers, nb_depots, vehicle_capacity, opening_route_cost, demands, capacity_depots,
            opening_depots_cost, distance_customers, distance_customers_depots)

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

# Computes the distance depot matrix
def compute_distance_depot(customers_x, customers_y, depot_x, depot_y, are_cost_double):
    nb_customers = len(customers_x)
    nb_depots = len(depot_x)
    distance_customers_depots = [[None for _ in range(nb_depots)] for _ in range(nb_customers)]
    for i in range(nb_customers):
        for d in range(nb_depots):
            dist = compute_dist(customers_x[i], depot_x[d], customers_y[i], depot_y[d], are_cost_double)
            distance_customers_depots[i][d] = dist
    return distance_customers_depots


def compute_dist(xi, xj, yi, yj, are_cost_double):
    dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
    if are_cost_double == 0:
        dist = math.ceil(100 * dist)
    return dist


def read_input_lrp(filename):
    if filename.endswith(".dat"):
        return read_input_lrp_dat(filename)
    else:
        raise Exception("Wrong file type. Please try again with a .dat file.")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python location_routing_problem.py input_file [output_file] [time_limit]")
        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"

    main(instance_file, str_time_limit, sol_file)
Compilation / Execution (Windows)
cl /EHsc location_routing_problem.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver110.lib
location_routing_problem instances/coordChrist100.dat
Compilation / Execution (Linux)
g++ location_routing_problem.cpp -I/opt/localsolver_11_0/include -llocalsolver110 -lpthread -o location_routing_problem
./location_routing_problem instances/coordChrist100.dat
/*************location_routing_problem.cpp*************/

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

using namespace localsolver;
using namespace std;

class LocationRoutingProblem {
  public:
    // LocalSolver
    LocalSolver localsolver;

    // Number of customers
    int nbCustomers;

    // Customers coordinates
    vector<int> xCustomers;
    vector<int> yCustomers;

    // Customers demands
    vector<double> demands;

    // Number of depots
    int nbDepots;

    // Depots coordinates
    vector<double> xDepots;
    vector<double> yDepots;

    // Capacity of depots
    vector<double> depotsCapacity;

    // Cost of opening a depot
    vector<double> openingDepotsCost;

    // Number of trucks
    int nbTrucks;

    // Capacity of trucks
    int truckCapacity;

    // Cost of opening a route
    int openingRouteCost;

    // Is the route used ?
    vector<LSExpression> sequenceUsed;

    // What is the depot of the route ?
    vector<LSExpression> associatedDepot;

    // Distance matrixes
    vector<vector<double>> distCustomers;
    vector<vector<double>> distCustomersDepots;

    int areCostDouble;

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

    // Sum of all the costs
    LSExpression totalCost;

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

    void solve(const char* limit) {
        // Declares the optimization model.
        LSModel m = localsolver.getModel();
        int minNbTrucks = ceil(accumulate(demands.begin(), demands.end(), 0) / truckCapacity);
        nbTrucks = ceil(1.5 * minNbTrucks);

        // A sequence is represented as a list containing the customers in the order they are visited
        customersSequences.resize(nbTrucks);
        for (int i = 0; i < nbTrucks; i++) {
            customersSequences[i] = m.listVar(nbCustomers);
        }
        // All customers should be assigned to a sequence
        m.constraint(m.partition(customersSequences.begin(), customersSequences.end()));

        // A depot is represented as a set containing the associated customersSequences
        depots.resize(nbDepots);
        for (int d = 0; d < nbDepots; d++) {
            depots[d] = m.setVar(nbTrucks);
        }
        // All the customersSequences should be assigned to a depot
        m.constraint(m.partition(depots.begin(), depots.end()));

        LSExpression quantityServed = m.array();
        LSExpression demandsArray = m.array(demands.begin(), demands.end());
        vector<LSExpression> sequenceLength;
        vector<LSExpression> sequenceCost;

        sequenceLength.resize(nbTrucks);
        sequenceUsed.resize(nbTrucks);
        sequenceCost.resize(nbTrucks);
        associatedDepot.resize(nbTrucks);
        LSExpression distCustomersArray = m.array();
        LSExpression distCustomersDepotsArray = m.array();

        // Create distance as an array to be able to acces it with an "at" operator
        for (int i = 0; i < nbCustomers; i++) {
            distCustomersArray.addOperand(m.array(distCustomers[i].begin(), distCustomers[i].end()));
            distCustomersDepotsArray.addOperand(m.array(distCustomersDepots[i].begin(), distCustomersDepots[i].end()));
        }

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

            // A sequence is used if it serves at least one customer
            sequenceUsed[r] = c > 0;
            // The "find" function gets the depot assigned to the sequence
            associatedDepot[r] = m.find(m.array(depots.begin(), depots.end()), r);

            LSExpression demandSelector = m.lambdaFunction([&](LSExpression i) { return demandsArray[i]; });
            quantityServed.addOperand(m.sum(sequence, demandSelector));
            // The quantity needed in each sequence must not exceed the vehicle capacity
            m.constraint(quantityServed[r] <= truckCapacity);

            LSExpression distSelector = m.lambdaFunction(
                [&](LSExpression i) { return m.at(distCustomersArray, sequence[i], sequence[i + 1]); });

            sequenceLength[r] = m.iif(sequenceUsed[r],
                                      m.at(distCustomersDepotsArray, sequence[0], associatedDepot[r]) +
                                          m.at(distCustomersDepotsArray, sequence[c - 1], associatedDepot[r]),
                                      0) +
                                m.sum(m.range(0, c - 1), distSelector);

            // The sequence cost is the sum of the opening cost and the sequence length
            sequenceCost[r] = sequenceUsed[r] * openingRouteCost + sequenceLength[r];
        }

        vector<LSExpression> depotCost;
        depotCost.resize(nbDepots);
        for (int d = 0; d < nbDepots; d++) {
            // A depot is open if at least a sequence starts from there
            depotCost[d] = openingDepotsCost[d] * (m.count(depots[d]) > 0);

            LSExpression depotSelect = m.lambdaFunction([&](LSExpression r) { return quantityServed[r]; });
            LSExpression depotQuantity = m.sum(depots[d], depotSelect);
            // The total demand served by a depot must not exceed its capacity
            m.constraint(depotQuantity <= depotsCapacity[d]);
        }
        LSExpression depotsCost = m.sum(depotCost.begin(), depotCost.end());
        LSExpression routingCost = m.sum(sequenceCost.begin(), sequenceCost.end());
        totalCost = routingCost + depotsCost;

        m.minimize(totalCost);
        m.close();

        localsolver.getParam().setTimeLimit(atoi(limit));
        localsolver.solve();
    }

    void writeSolution(const char* inFile, const string& solFile) {
        ofstream file;
        file.exceptions(ofstream::failbit | ofstream::badbit);
        file.open(solFile.c_str());
        file << "File name: " << inFile << "; total cost = " << totalCost.getDoubleValue() << endl;
        for (int r = 0; r < nbTrucks; r++) {
            if (sequenceUsed[r].getValue()) {
                file << "Sequence " << r << ", assigned to depot " << associatedDepot[r].getValue() << " : ";
                LSCollection customersCollection = customersSequences[r].getCollectionValue();
                for (lsint i = 0; i < customersCollection.count(); i++) {
                    file << customersCollection[i] << " ";
                }
                file << endl;
            }
        }
    }

  private:
    void readInputLrp(const char* fileName) {
        string file = fileName;
        ifstream infile(file.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }
        infile >> nbCustomers;
        xCustomers.resize(nbCustomers);
        yCustomers.resize(nbCustomers);
        demands.resize(nbCustomers);
        distCustomers.resize(nbCustomers);
        distCustomersDepots.resize(nbCustomers);
        infile >> nbDepots;
        xDepots.resize(nbDepots);
        yDepots.resize(nbDepots);
        depotsCapacity.resize(nbDepots);
        openingDepotsCost.resize(nbDepots);
        for (int i = 0; i < nbDepots; i++) {
            infile >> xDepots[i];
            infile >> yDepots[i];
        }
        for (int i = 0; i < nbCustomers; i++) {
            infile >> xCustomers[i];
            infile >> yCustomers[i];
        }
        infile >> truckCapacity;
        for (int i = 0; i < nbDepots; i++) {
            infile >> depotsCapacity[i];
        }
        for (int i = 0; i < nbCustomers; i++) {
            infile >> demands[i];
        }
        vector<double> tempOpeningCostDepots;
        tempOpeningCostDepots.resize(nbDepots);
        for (int i = 0; i < nbDepots; i++) {
            infile >> tempOpeningCostDepots[i];
        }
        int tempOpeningCostRoute;
        infile >> tempOpeningCostRoute;
        infile >> areCostDouble;
        infile.close();
        if (areCostDouble == 1) {
            openingRouteCost = tempOpeningCostRoute;
            for (int i = 0; i < nbDepots; i++) {
                openingDepotsCost[i] = tempOpeningCostDepots[i];
            }
        } else {
            openingRouteCost = round(tempOpeningCostRoute);
            for (int i = 0; i < nbDepots; i++) {
                openingDepotsCost[i] = round(tempOpeningCostDepots[i]);
            }
        }
        computeDistanceMatrix();
    }

    void computeDistanceMatrix() {
        for (int i = 0; i < nbCustomers; i++) {
            distCustomers[i].resize(nbCustomers);
            distCustomersDepots[i].resize(nbDepots);
            for (int j = 0; j < nbCustomers; j++) {
                distCustomers[i][j] =
                    computeDist(xCustomers[i], yCustomers[i], xCustomers[j], yCustomers[j], areCostDouble);
            }
            for (int d = 0; d < nbDepots; d++) {
                distCustomersDepots[i][d] =
                    computeDist(xCustomers[i], yCustomers[i], xDepots[d], yDepots[d], areCostDouble);
            }
        }
    }

    double computeDist(int xi, int yi, int xj, int yj, int areCostDouble) {
        double dist = sqrt(pow(xi - xj, 2) + pow(yi - yj, 2));
        if (areCostDouble == 0) {
            dist = ceil(dist * 100);
        }
        return dist;
    }
};

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

    try {
        LocationRoutingProblem model;
        model.readInstance(instanceFile);
        model.solve(strTimeLimit);
        if (solFile != NULL)
            model.writeSolution(instanceFile, solFile);
    } catch (const std::exception& e) {
        std::cerr << "An error occured: " << e.what() << endl;
    }

    return 0;
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc LocationRoutingProblem.cs /reference:localsolvernet.dll
LocationRoutingProblem instances\coordChrist100.dat
/************LocationRoutingProblem.cs*************/

using System;
using System.IO;
using System.Collections.Generic;
using localsolver;

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

    // customers data
    int nbCustomers;

    // Customers coordinates
    int[] xCustomers;
    int[] yCustomers;

    // Customers demands
    int[] demands;

    // Depots data
    int nbDepots;

    // Depots coordinates
    int[] xDepots;
    int[] yDepots;

    // Capacity of depots
    int[] depotsCapacity;
    double[] openingDepotsCost;

    // Number of trucks
    int nbTrucks;

    // Capacity of a truck
    int truckCapacity;

    // Cost of opening a route
    int openingRouteCost;

    // Is the route used ?
    LSExpression[] sequenceUsed;

    // What is the depot of the route ?
    LSExpression[] associatedDepot;

    // Distance matrices
    double[][] distCustomers;
    double[][] distCustomersDepots;

    int areCostDouble;

    // Decision variables
    LSExpression[] customersSequences;
    LSExpression[] depots;

    // Sum of all the costs
    LSExpression totalCost;

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

    // Reads instance data.
    void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] line;
            line = ReadNextLine(input);
            nbCustomers = int.Parse(line[0]);
            xCustomers = new int[nbCustomers];
            yCustomers = new int[nbCustomers];
            demands = new int[nbCustomers];
            distCustomers = new double[nbCustomers][];

            line = ReadNextLine(input);
            nbDepots = int.Parse(line[0]);
            xDepots = new int[nbDepots];
            yDepots = new int[nbDepots];
            depotsCapacity = new int[nbDepots];
            openingDepotsCost = new double[nbDepots];
            distCustomersDepots = new double[nbCustomers][];

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

            for (int i = 0; i < nbCustomers; i++)
            {
                line = ReadNextLine(input);
                xCustomers[i] = int.Parse(line[0]);
                yCustomers[i] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            truckCapacity = int.Parse(line[0]);
            line = ReadNextLine(input);

            for (int i = 0; i < nbDepots; i++)
            {
                line = ReadNextLine(input);
                depotsCapacity[i] = int.Parse(line[0]);
            }
            line = ReadNextLine(input);
            for (int i = 0; i < nbCustomers; i++)
            {
                line = ReadNextLine(input);
                demands[i] = int.Parse(line[0]);
            }
            line = ReadNextLine(input);

            double[] tempOpeningCostDepots = new double[nbDepots];
            for (int i = 0; i < nbDepots; i++)
            {
                line = ReadNextLine(input);
                tempOpeningCostDepots[i] = double.Parse(
                    line[0],
                    System.Globalization.CultureInfo.InvariantCulture
                );
            }
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            double tempOpeningCostRoute = double.Parse(line[0]);
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            areCostDouble = int.Parse(line[0]);

            if (areCostDouble == 1)
            {
                openingRouteCost = (int)tempOpeningCostRoute;
                for (int i = 0; i < nbDepots; i++)
                {
                    openingDepotsCost[i] = tempOpeningCostDepots[i];
                }
            }
            else
            {
                openingRouteCost = (int)Math.Round(tempOpeningCostRoute);
                for (int i = 0; i < nbDepots; i++)
                {
                    openingDepotsCost[i] = Math.Round(tempOpeningCostDepots[i]);
                }
            }
            DistanceMatrices();
        }
    }

    void DistanceMatrices()
    {
        for (int i = 0; i < nbCustomers; i++)
        {
            distCustomers[i] = new double[nbCustomers];
            for (int j = 0; j < nbCustomers; j++)
            {
                distCustomers[i][j] = ComputeDistance(
                    xCustomers[i],
                    yCustomers[i],
                    xCustomers[j],
                    yCustomers[j]
                );
            }
            distCustomersDepots[i] = new double[nbDepots];
            for (int d = 0; d < nbDepots; d++)
            {
                distCustomersDepots[i][d] = ComputeDistance(
                    xCustomers[i],
                    yCustomers[i],
                    xDepots[d],
                    yDepots[d]
                );
            }
        }
    }

    double ComputeDistance(int xi, int yi, int xj, int yj)
    {
        double dist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        if (areCostDouble == 0)
        {
            dist = Math.Ceiling(dist * 100);
        }
        return dist;
    }

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

    void Solve(int limit)
    {
        LSModel m = localsolver.GetModel();
        int totalDemand = 0;
        foreach (int dem in demands)
        {
            totalDemand += dem;
        }
        int minNbTrucks = (int)Math.Ceiling((double)totalDemand / truckCapacity);
        nbTrucks = (int)Math.Ceiling(1.5 * minNbTrucks);

        customersSequences = new LSExpression[nbTrucks];
        depots = new LSExpression[nbDepots];

        // A sequence is represented as a list containing the customers in the order they are visited
        for (int r = 0; r < nbTrucks; r++)
        {
            customersSequences[r] = m.List(nbCustomers);
        }
        // All customers should be assigned to a sequence
        m.Constraint(m.Partition(customersSequences));
        // A depot is represented as a set containing the associated sequences
        for (int d = 0; d < nbDepots; d++)
        {
            depots[d] = m.Set(nbTrucks);
        }
        // All the sequences should be assigned to a depot
        m.Constraint(m.Partition(depots));

        LSExpression demandsArray = m.Array(demands);
        LSExpression distCustomersArray = m.Array(distCustomers);
        LSExpression distCustomersDepotsArray = m.Array(distCustomersDepots);
        sequenceUsed = new LSExpression[nbTrucks];
        LSExpression[] sequenceCost = new LSExpression[nbTrucks];
        LSExpression[] sequenceLength = new LSExpression[nbTrucks];
        associatedDepot = new LSExpression[nbTrucks];
        LSExpression quantityServed = m.Array();

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

            // A sequence is used if it serves at least one customer
            sequenceUsed[r] = m.Gt(c, 0);
            // The "find" function gets the depot that is assigned to the sequence
            associatedDepot[r] = m.Find(m.Array(depots), r);

            LSExpression demandSelector = m.LambdaFunction(i => demandsArray[i]);
            quantityServed.AddOperand(m.Sum(sequence, demandSelector));
            // The quantity needed in each sequence must not exceed the vehicle capacity
            m.Constraint(quantityServed[r] <= truckCapacity);

            LSExpression distSelector = m.LambdaFunction(
                i => m.At(distCustomersArray, m.At(sequence, i), m.At(sequence, i + 1))
            );
            
            sequenceLength[r] =
                m.Sum(m.Range(0, c - 1), distSelector)
                + m.If(
                    sequenceUsed[r],
                    m.At(distCustomersDepotsArray, m.At(sequence, 0), associatedDepot[r])
                        + m.At(distCustomersDepotsArray, m.At(sequence, c - 1), associatedDepot[r]),
                    0
                );
            // The sequence cost is the sum of the opening cost and the sequence length
            sequenceCost[r] = sequenceLength[r] + openingRouteCost * sequenceUsed[r];
        }
        LSExpression[] depotCost = new LSExpression[nbDepots];
        for (int d = 0; d < nbDepots; d++)
        {
            // A depot is open if at least a sequence starts from there
            depotCost[d] = openingDepotsCost[d] * (m.Count(depots[d]) > 0);
            LSExpression depotSelect = m.LambdaFunction(r => m.At(quantityServed, r));
            LSExpression depotQuantity = m.Sum(depots[d], depotSelect);
            // The total demand served by a depot must not exceed its capacity
            m.Constraint(depotQuantity <= depotsCapacity[d]);
        }
        LSExpression depotsCost = m.Sum(depotCost);
        LSExpression routingCost = m.Sum(sequenceCost);
        totalCost = routingCost + depotsCost;

        m.Minimize(totalCost);
        m.Close();

        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.Solve();
    }

    void WriteSolution(string infile, string outfile)
    {
        using (StreamWriter output = new StreamWriter(outfile))
        {
            output.WriteLine(
                "File name: " + infile + "; total cost = " + totalCost.GetDoubleValue()
            );
            for (int r = 0; r < nbTrucks; r++)
            {
                if (sequenceUsed[r].GetValue() != 0)
                {
                    output.Write(
                        "Route " + r + ", assigned to depot " + associatedDepot[r].GetValue() + ": "
                    );
                    LSCollection customersCollection = customersSequences[r].GetCollectionValue();
                    for (int i = 0; i < customersCollection.Count(); i++)
                    {
                        output.Write(customersCollection[i] + " ");
                    }
                    output.WriteLine();
                }
            }
        }
    }

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

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

        using (LocationRoutingProblem model = new LocationRoutingProblem())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (strOutput != null)
                model.WriteSolution(instanceFile, strOutput);
        }
    }
}
Compilation / Execution (Windows)
javac LocationRoutingProblem.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. LocationRoutingProblem instances\coordChrist100.dat
Compilation / Execution (Linux)
javac LocationRoutingProblem.java -cp /opt/localsolver_11_0/bin/localsolver.jar
java -cp /opt/localsolver_11_0/bin/localsolver.jar:. LocationRoutingProblem instances/coordChrist100.dat
/*************LocationRoutingProblem.java*************/

import localsolver.*;
import java.io.*;
import java.io.File; // Import the File class
import java.io.FileNotFoundException; // Import this class to handle errors
import java.util.Scanner; // Import the Scanner class to read text files

public class LocationRoutingProblem {
    // Solver
    private final LocalSolver localsolver;

    // Number of customers
    private int nbCustomers;

    // Customers coordinates
    private int[] xCustomers;
    private int[] yCustomers;

    // Customers demands
    private int[] demands;

    // Depots data
    private int nbDepots;

    // Depots coordinates
    private int[] xDepots;
    private int[] yDepots;

    // Capacity of depots
    private int[] depotsCapacity;

    // Cost of opening a depot
    private double[] openingDepotsCost;

    // Number of trucks
    private int nbTrucks;

    // Capacity of trucks
    private int truckCapacity;

    // Cost of opening a route
    private int openingRouteCost;

    // Is the sequence used ?
    private LSExpression[] sequenceUsed;

    // What is the depot of the sequence ?
    private LSExpression[] associatedDepot;

    // Distance matrixes
    private double[][] distCustomers;
    private double[][] distCustomersDepots;

    private int areCostDouble;

    // Decision variables
    private LSExpression[] customersSequences;
    private LSExpression[] depots;

    // Objective value
    private LSExpression totalCost;

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

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

        int totalDemand = 0;
        for (int d : demands) {
            totalDemand += d;
        }
        int minNbTrucks = (int) Math.ceil(totalDemand / truckCapacity);
        int nbTrucks = (int) Math.ceil(1.5 * minNbTrucks);

        customersSequences = new LSExpression[nbTrucks];
        depots = new LSExpression[nbDepots];

        // A sequence is represented as a list containing the customers in the order
        // they are visited
        for (int i = 0; i < nbTrucks; i++) {
            customersSequences[i] = m.listVar(nbCustomers);
        }
        // All customers should be assigned to a sequence
        m.constraint(m.partition(customersSequences));
        // A depot is represented as a set containing the associated customersSequences
        for (int d = 0; d < nbDepots; d++) {
            depots[d] = m.setVar(nbTrucks);
        }
        // All the customersSequences should be assigned to a depot
        m.constraint(m.partition(depots));

        LSExpression demandsArray = m.array(demands);
        LSExpression distCustomersArray = m.array(distCustomers);
        LSExpression distCustomersDepotsArray = m.array(distCustomersDepots);
        sequenceUsed = new LSExpression[nbTrucks];
        LSExpression[] sequenceCost = new LSExpression[nbTrucks];
        LSExpression[] sequenceLength = new LSExpression[nbTrucks];
        associatedDepot = new LSExpression[nbTrucks];
        LSExpression quantityServed = m.array();

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

            // A sequence is used if it serves at least one customer
            sequenceUsed[r] = m.gt(c, 0);
            // The "find" function gets the depot that is assigned to the sequence
            associatedDepot[r] = m.find(m.array(depots), r);

            LSExpression demandSelector = m.lambdaFunction(i -> m.at(demandsArray, i));
            quantityServed.addOperand(m.sum(sequence, demandSelector));
            // The quantity needed in each sequence must not exceed the vehicle capacity
            m.constraint(m.leq(m.at(quantityServed, r), truckCapacity));

            LSExpression distSelector = m
                    .lambdaFunction(i -> m.at(distCustomersArray, m.at(sequence, i), m.at(sequence, m.sum(i, 1))));
            
            sequenceLength[r] = m.sum(m.sum(m.range(0, m.sub(c, 1)), distSelector),
                    m.iif(sequenceUsed[r],
                            m.sum(m.at(distCustomersDepotsArray, m.at(sequence, 0), associatedDepot[r]), 
                                  m.at(distCustomersDepotsArray, m.at(sequence, m.sub(c, 1)), associatedDepot[r])),
                            0));
            // The sequence cost is the sum of the opening cost and the sequence length
            sequenceCost[r] = m.sum(sequenceLength[r], m.prod(openingRouteCost, sequenceUsed[r]));
        }

        LSExpression[] depotCost = new LSExpression[nbDepots];
        for (int d = 0; d < nbDepots; d++) {
            // A depot is open if at least a sequence starts from there
            depotCost[d] = m.prod(openingDepotsCost[d], m.gt(m.count(depots[d]), 0));

            LSExpression depotSelect = m.lambdaFunction(r -> m.at(quantityServed, r));
            LSExpression depotQuantity = m.sum(depots[d], depotSelect);
            // The total demand served by a depot must not exceed its capacity
            m.constraint(m.leq(depotQuantity, depotsCapacity[d]));
        }
        LSExpression depotsCost = m.sum(depotCost);
        LSExpression routingCost = m.sum(sequenceCost);
        totalCost = m.sum(routingCost, depotsCost);

        m.minimize(totalCost);
        m.close();

        localsolver.getParam().setTimeLimit(limit);
        localsolver.solve();
    }

    private void writeSolution(String infile, String outfile) throws IOException {
        try (PrintWriter output = new PrintWriter(outfile)) {
            output.println("File name: " + infile + "; total cost = " + totalCost.getDoubleValue());
            for (int r = 0; r < nbTrucks; r++) {
                if (sequenceUsed[r].getIntValue() != 0) {
                    output.print("Route " + r + ", assigned to depot " + associatedDepot[r].getIntValue() + " : ");
                    LSCollection customersCollection = customersSequences[r].getCollectionValue();
                    for (int i = 0; i < customersCollection.count(); i++) {
                        output.print(customersCollection.get(i) + " ");
                    }
                    output.println();
                }
            }
        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    private void readInstanceLrp(String fileName) throws IOException {
        try (Scanner infile = new Scanner(new File(fileName))) {
            nbCustomers = infile.nextInt();
            xCustomers = new int[nbCustomers];
            yCustomers = new int[nbCustomers];
            demands = new int[nbCustomers];
            distCustomers = new double[nbCustomers][];
            distCustomersDepots = new double[nbCustomers][];
            nbDepots = infile.nextInt();
            xDepots = new int[nbDepots];
            yDepots = new int[nbDepots];
            depotsCapacity = new int[nbDepots];
            openingDepotsCost = new double[nbDepots];

            for (int i = 0; i < nbDepots; i++) {
                xDepots[i] = infile.nextInt();
                yDepots[i] = infile.nextInt();
            }
            for (int i = 0; i < nbCustomers; i++) {
                xCustomers[i] = infile.nextInt();
                yCustomers[i] = infile.nextInt();
            }
            truckCapacity = infile.nextInt();
            for (int i = 0; i < nbDepots; i++) {
                depotsCapacity[i] = infile.nextInt();
            }
            for (int i = 0; i < nbCustomers; i++) {
                demands[i] = infile.nextInt();
            }
            double[] tempOpeningCostDepots = new double[nbDepots];
            for (int i = 0; i < nbDepots; i++) {
                if (infile.hasNext()) {
                    tempOpeningCostDepots[i] = Double.parseDouble(infile.next());
                } else if (infile.hasNextInt()) {
                    tempOpeningCostDepots[i] = infile.nextInt();
                }
            }
            int tempOpeningCostRoute = infile.nextInt();
            areCostDouble = infile.nextInt();
            if (areCostDouble == 1) {
                openingRouteCost = tempOpeningCostRoute;
                for (int i = 0; i < tempOpeningCostDepots.length; i++) {
                    openingDepotsCost[i] = tempOpeningCostDepots[i];
                }
            } else {
                openingRouteCost = Math.round(tempOpeningCostRoute);
                for (int i = 0; i < tempOpeningCostDepots.length; i++) {
                    openingDepotsCost[i] = Math.round(tempOpeningCostDepots[i]);
                }
            }
            computeDistanceMatrix();
        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    void computeDistanceMatrix() {
        for (int i = 0; i < nbCustomers; i++) {
            distCustomers[i] = new double[nbCustomers];
            for (int j = 0; j < nbCustomers; j++) {
                distCustomers[i][j] = computeDist(xCustomers[i], xCustomers[j], yCustomers[i], yCustomers[j],
                        areCostDouble);
            }
            distCustomersDepots[i] = new double[nbDepots];
            for (int d = 0; d < nbDepots; d++) {
                distCustomersDepots[i][d] = computeDist(xCustomers[i], xDepots[d], yCustomers[i], yDepots[d],
                        areCostDouble);
            }
        }
    }

    private double computeDist(int xi, int xj, int yi, int yj, int areCostDouble) {
        double dist = Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
        if (areCostDouble == 0) {
            dist = Math.ceil(100 * dist);
        }
        return dist;
    }

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

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

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