View all examples

Location Routing Problem (LRP)

Problem

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, several depot locations are available, each of them having its own opening cost and capacity. The vehicles start and end their routes at the same depot. Each customer must be served by exactly one vehicle, and the total demand served by each vehicle and each depot must not exceed its capacity. The objective is to minimize the total cost, which is the sum of the route lengths and opening costs for the depots and the routes.

Principles learned

Data

The instances we provide 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
  • A Boolean value, indicating whether costs should be rounded.

Program

The Hexaly model for the Location Routing Problem (LRP) uses list decision variables. For each truck, we define a list variable representing the sequence of customers it visits. Using a ‘partition’ constraint on all the lists, we ensure that each customer is served by exactly one truck. In addition, we also use set decision variables. For each depot, we define a set variable representing the trucks starting and ending their routes at this depot. Using another ‘partition’ constraint on the sets, we ensure that each truck is associated with a depot.

The total quantity delivered by each truck is computed with a lambda function to apply the ‘sum’ operator over all visited customers. Note that the number of terms in this sum varies during the search, along with the size of the list. We can then constrain this quantity to be lower than the truck capacity.

Using the ‘find‘ operator, we can then retrieve the index of the depot associated with each truck. This allows us to compute the distance from each truck’s depot to its first customer, and from its last customer back to its depot. To compute the total length of each route, we also need to know the distance from one customer to the next: we use another lambda function to sum the distances along the route.

A route is considered open if it serves at least one customer, and a depot is considered open if at least one non-empty sequence uses it. These conditions are computed thanks to the ‘count’ operator.

Finally, we can compute the objective function, by summing the total length of the routes and the opening costs of each open route and depot.

Download
Execution
hexaly location_routing_problem.hxm inFileName=instances/coordChrist100.dat [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data */
function input() {

    local usage = "Usage: hexaly location_routing_problem.hxm "
            + "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";
    if (inFileName == nil) throw usage;

    readInputLrp();
    computeDistances();

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

/* Declare the optimization model */
function model() {
    // A route is represented as a list containing the customers in the order they are
    // visited
    customersSequences[r in 0...nbTrucks] <- list(nbCustomers);
    // A depot is represented as a set containing the associated customersSequences
    depots[d in 0...nbDepots] <- 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] {
        local sequence <- customersSequences[r];
        local c <- count(sequence);

        quantityServed[r] <- sum(sequence, j => demands[j]);
        // 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-1, 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] {
        // 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](isDepotOpen[d] * openingCostDepots[d]);
    routingCost <- sum[r in 0...nbTrucks](routeCost[r]);

    totalCost <- routingCost + depotsCost;
    minimize totalCost;
}

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

/* Write the solution in a file */
function output() {
    if (solFileName == nil) return;
    local resFile = io.openWrite(solFileName);
    resFile.println("File name: ", inFileName + "; totalCost = " + totalCost.value);
    for [r in 0...nbTrucks] {
        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[d in 0...nbDepots][l in 0...2] = inFile.readDouble();
    coordinatesCustomers[c in 0...nbCustomers][l in 0...2] = inFile.readDouble();

    truckCapacity = inFile.readInt();
    depotsCapacity[d in 0...nbDepots] = inFile.readInt();

    demands[c in 0...nbCustomers] = inFile.readInt();

    local tempOpeningCostDepots[d in 0...nbDepots] = inFile.readDouble();
    local tempOpeningCostRoute = inFile.readDouble();

    areCostsDouble = inFile.readInt();

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

function computeDistances() {
    local tempDistanceCustomers[c0 in 0...nbCustomers][c1 in 0...nbCustomers] =
            sqrt(pow((coordinatesCustomers[c0][0] - coordinatesCustomers[c1][0]), 2)
            + pow((coordinatesCustomers[c0][1] - coordinatesCustomers[c1][1]), 2));
    local tempDistanceCustomersDepots[c in 0...nbCustomers][d in 0...nbDepots] =
            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][c1 in 0...nbCustomers] =
                tempDistanceCustomers[c0][c1];
        distCustomersDepots[c in 0...nbCustomers][d in 0...nbDepots] =
                tempDistanceCustomersDepots[c][d];
    } else {
        distCustomers[c0 in 0...nbCustomers][c1 in 0...nbCustomers] =
                ceil(100 * tempDistanceCustomers[c0][c1]);
        distCustomersDepots[c in 0...nbCustomers][d in 0...nbDepots] =
                ceil(100 * tempDistanceCustomersDepots[c][d]);
    }
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python location_routing_problem.py instances\coordChrist100.da
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python location_routing_problem.py instances/coordChrist100.dat
import hexaly.optimizer
import sys
import math


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


def main(instance_file, str_time_limit, sol_file):
    #
    # Read instance data
    #
    nb_customers, nb_depots, vehicle_capacity, opening_route_cost, demands_data, \
        capacity_depots, opening_depots_cost, dist_matrix_data, dist_depots_data = \
        read_input_lrp(instance_file)

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

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

        # A route 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 assigned to a route
        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 assigned to a depot
        m.constraint(m.partition(depots))

        route_costs = [None] * nb_trucks
        sequence_used = [None] * nb_trucks
        dist_routes = [None] * nb_trucks
        associated_depot = [None] * nb_trucks

        # Create Hexaly arrays to be able to access them with "at" operators
        demands = m.array(demands_data)
        dist_matrix = m.array()
        dist_depots = m.array()
        quantity_served = m.array()
        for i in range(nb_customers):
            dist_matrix.add_operand(m.array(dist_matrix_data[i]))
            dist_depots.add_operand(m.array(dist_depots_data[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_lambda = m.lambda_function(lambda j: demands[j])
            quantity_served.add_operand(m.sum(sequence, demand_lambda))
            m.constraint(quantity_served[r] <= vehicle_capacity)

            # Distance traveled by each truck
            dist_lambda = m.lambda_function(
                lambda i: m.at(dist_matrix, sequence[i], sequence[i + 1]))
            depot = associated_depot[r]
            dist_routes[r] = m.sum(m.range(0, c - 1), dist_lambda) + m.iif(
                sequence_used[r],
                m.at(dist_depots, sequence[0], depot)
                + m.at(dist_depots, sequence[c - 1], depot),
                0)

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

        depot_cost = [None] * nb_depots
        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_lambda = m.lambda_function(lambda r: quantity_served[r])
            depot_quantity = m.sum(depots[d], depot_lambda)
            m.constraint(depot_quantity <= capacity_depots[d])

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

        m.minimize(totalCost)

        m.close()

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

        optimizer.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

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

# Compute 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("Unknown file format")


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%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
location_routing_problem instances\coordChrist100.dat
Compilation / Execution (Linux)
g++ location_routing_problem.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o location_routing_problem
./location_routing_problem instances/coordChrist100.dat
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <numeric>
#include <vector>

using namespace hexaly;
using namespace std;

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

    // Number of customers
    int nbCustomers;

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

    // Customers demands
    vector<double> demandsData;

    // 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<HxExpression> sequenceUsed;

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

    // Distance matrixes
    vector<vector<double>> distMatrixData;
    vector<vector<double>> distDepotsData;

    int areCostDouble;

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

    // Sum of all the costs
    HxExpression totalCost;

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

    void solve(const char* limit) {
        // Declare the optimization model
        HxModel m = optimizer.getModel();
        int minNbTrucks = ceil(accumulate(demandsData.begin(), demandsData.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()));

        vector<HxExpression> distRoutes;
        vector<HxExpression> routeCosts;
        distRoutes.resize(nbTrucks);
        sequenceUsed.resize(nbTrucks);
        routeCosts.resize(nbTrucks);
        associatedDepot.resize(nbTrucks);

        // Create Hexaly arrays to be able to access them with "at" operators
        HxExpression quantityServed = m.array();
        HxExpression demands = m.array(demandsData.begin(), demandsData.end());
        HxExpression distMatrix = m.array();
        HxExpression distDepots = m.array();
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrix.addOperand(m.array(distMatrixData[i].begin(), distMatrixData[i].end()));
            distDepots.addOperand(m.array(distDepotsData[i].begin(), distDepotsData[i].end()));
        }

        for (int r = 0; r < nbTrucks; ++r) {
            HxExpression sequence = customersSequences[r];
            HxExpression 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);

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

            HxExpression distLambda =
                m.lambdaFunction([&](HxExpression i) { return m.at(distMatrix, sequence[i], sequence[i + 1]); });

            distRoutes[r] = m.iif(sequenceUsed[r],
                                  m.at(distDepots, sequence[0], associatedDepot[r]) +
                                      m.at(distDepots, sequence[c - 1], associatedDepot[r]),
                                  0) +
                            m.sum(m.range(0, c - 1), distLambda);

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

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

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

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

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

    /* Write the solution in a file */
    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() << " : ";
                HxCollection customersCollection = customersSequences[r].getCollectionValue();
                for (hxint 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);
        demandsData.resize(nbCustomers);
        distMatrixData.resize(nbCustomers);
        distDepotsData.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 >> demandsData[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) {
            distMatrixData[i].resize(nbCustomers);
            distDepotsData[i].resize(nbDepots);
            for (int j = 0; j < nbCustomers; ++j) {
                distMatrixData[i][j] =
                    computeDist(xCustomers[i], yCustomers[i], xCustomers[j], yCustomers[j], areCostDouble);
            }
            for (int d = 0; d < nbDepots; ++d) {
                distDepotsData[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 %HX_HOME%\bin\Hexaly.NET.dll .
csc LocationRoutingProblem.cs /reference:Hexaly.NET.dll
LocationRoutingProblem instances\coordChrist100.dat
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

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

    // Customers demands
    int[] demandsData;

    // Number of depots
    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 ?
    HxExpression[] sequenceUsed;

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

    // Distance matrixes
    double[][] distMatrixData;
    double[][] distDepotsData;

    int areCostDouble;

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

    // Sum of all the costs
    HxExpression totalCost;

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

    /* Read 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];
            demandsData = new int[nbCustomers];
            distMatrixData = 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];
            distDepotsData = 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);
                demandsData[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]);
            }
            DistanceMatrixes();
        }
    }

    void DistanceMatrixes()
    {
        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i] = new double[nbCustomers];
            for (int j = 0; j < nbCustomers; ++j)
            {
                distMatrixData[i][j] = ComputeDistance(
                    xCustomers[i],
                    yCustomers[i],
                    xCustomers[j],
                    yCustomers[j]
                );
            }
            distDepotsData[i] = new double[nbDepots];
            for (int d = 0; d < nbDepots; ++d)
            {
                distDepotsData[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 (optimizer != null)
            optimizer.Dispose();
    }

    void Solve(int limit)
    {
        // Declare the optimization model
        HxModel m = optimizer.GetModel();
        int totalDemand = 0;
        foreach (int dem in demandsData)
            totalDemand += dem;
        int minNbTrucks = (int)Math.Ceiling((double)totalDemand / truckCapacity);
        nbTrucks = (int)Math.Ceiling(1.5 * minNbTrucks);

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

        // A route 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 route
        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));

        // Create HexalyOptimizer arrays to be able to access them with "at" operators
        HxExpression demands = m.Array(demandsData);
        HxExpression distMatrix = m.Array(distMatrixData);
        HxExpression distDepots = m.Array(distDepotsData);
        sequenceUsed = new HxExpression[nbTrucks];
        HxExpression[] routeCosts = new HxExpression[nbTrucks];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        associatedDepot = new HxExpression[nbTrucks];
        HxExpression quantityServed = m.Array();

        for (int r = 0; r < nbTrucks; ++r)
        {
            HxExpression sequence = customersSequences[r];
            HxExpression 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);

            HxExpression demandLambda = m.LambdaFunction(j => demands[j]);
            quantityServed.AddOperand(m.Sum(sequence, demandLambda));
            // The quantity needed in each sequence must not exceed the vehicle capacity
            m.Constraint(quantityServed[r] <= truckCapacity);

            HxExpression distLambda = m.LambdaFunction(
                i => m.At(distMatrix, m.At(sequence, i), m.At(sequence, i + 1))
            );

            distRoutes[r] =
                m.Sum(m.Range(0, c - 1), distLambda)
                + m.If(
                    sequenceUsed[r],
                    m.At(distDepots, m.At(sequence, 0), associatedDepot[r])
                        + m.At(distDepots, m.At(sequence, c - 1), associatedDepot[r]),
                    0
                );
            // The sequence cost is the sum of the opening cost and the sequence length
            routeCosts[r] = distRoutes[r] + openingRouteCost * sequenceUsed[r];
        }
        HxExpression[] depotCost = new HxExpression[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);
            HxExpression depotLambda = m.LambdaFunction(r => m.At(quantityServed, r));
            HxExpression depotQuantity = m.Sum(depots[d], depotLambda);
            // The total demand served by a depot must not exceed its capacity
            m.Constraint(depotQuantity <= depotsCapacity[d]);
        }
        HxExpression depotsCost = m.Sum(depotCost);
        HxExpression routingCost = m.Sum(routeCosts);
        totalCost = routingCost + depotsCost;

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

        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    /* Write the solution in a file */
    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() + ": "
                    );
                    HxCollection 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 %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. LocationRoutingProblem instances\coordChrist100.dat
Compilation / Execution (Linux)
javac LocationRoutingProblem.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. LocationRoutingProblem instances/coordChrist100.dat
import com.hexaly.optimizer.*;
import java.io.*;
import java.io.File;
import java.io.FileNotFoundException;
import java.util.Scanner;

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

    // Number of customers
    private int nbCustomers;

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

    // Customers demands
    private int[] demandsData;

    // Number of depots
    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 HxExpression[] sequenceUsed;

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

    // Distance matrixes
    private double[][] distMatrixData;
    private double[][] distDepotsData;

    private int areCostDouble;

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

    // Objective value
    private HxExpression totalCost;

    private LocationRoutingProblem(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

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

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

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

        // A route 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 route
        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));

        HxExpression demands = m.array(demandsData);
        HxExpression distMatrix = m.array(distMatrixData);
        HxExpression distDepots = m.array(distDepotsData);
        sequenceUsed = new HxExpression[nbTrucks];
        HxExpression[] routeCosts = new HxExpression[nbTrucks];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        associatedDepot = new HxExpression[nbTrucks];
        HxExpression quantityServed = m.array();

        for (int r = 0; r < nbTrucks; ++r) {
            HxExpression sequence = customersSequences[r];
            HxExpression 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);

            HxExpression demandLambda = m.lambdaFunction(j -> m.at(demands, j));
            quantityServed.addOperand(m.sum(sequence, demandLambda));
            // The quantity needed in each sequence must not exceed the vehicle capacity
            m.constraint(m.leq(m.at(quantityServed, r), truckCapacity));

            HxExpression distLambda = m
                .lambdaFunction(i -> m.at(distMatrix, m.at(sequence, i), m.at(sequence, m.sum(i, 1))));

            distRoutes[r] = m.sum(m.sum(m.range(0, m.sub(c, 1)), distLambda),
                m.iif(sequenceUsed[r], m.sum(m.at(distDepots, m.at(sequence, 0), associatedDepot[r]),
                    m.at(distDepots, m.at(sequence, m.sub(c, 1)), associatedDepot[r])), 0));
            // The sequence cost is the sum of the opening cost and the sequence length
            routeCosts[r] = m.sum(distRoutes[r], m.prod(openingRouteCost, sequenceUsed[r]));
        }

        HxExpression[] depotCost = new HxExpression[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));

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

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

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

    /* Write the solution in a file */
    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() + " : ");
                    HxCollection 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];
            demandsData = new int[nbCustomers];
            distMatrixData = new double[nbCustomers][];
            distDepotsData = 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) {
                demandsData[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) {
            distMatrixData[i] = new double[nbCustomers];
            for (int j = 0; j < nbCustomers; ++j) {
                distMatrixData[i][j] = computeDist(xCustomers[i], xCustomers[j], yCustomers[i], yCustomers[j],
                    areCostDouble);
            }
            distDepotsData[i] = new double[nbDepots];
            for (int d = 0; d < nbDepots; ++d) {
                distDepotsData[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 (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            String instanceFile = args[0];
            String strOutfile = args.length > 1 ? args[1] : null;
            String strTimeLimit = args.length > 2 ? args[2] : "20";

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

Similar examples

Discover the ease of use and performance of Hexaly