Multi Depot Vehicle Routing Problem (MDVRP)

Problem

In the Multi Depot Vehicle Routing Problem (MDVRP), a fleet of delivery vehicles with uniform capacity must service customers with known demand for a single commodity and service time. There are several depot locations, and each truck starts and ends its route at the same depot. The trucks available at each depot share a common capacity and maximum route duration. Each customer must be served by exactly one vehicle, and the total demand and route duration for each vehicle must not exceed its capacity and maximum duration. The objective is to minimize the total traveled distance.

Principles learned

Add list decision variables to model the trucks’ sequences of customers
Define lambda functions to compute the traveled distance and service time for each truck

Data

The instances we provide come from the Cordeau_2011 instances. The format of the data files is as follows:

First line: data type (number 2 for the MDVRP), number of trucks per depot, number of customers, number of depots
For each depot: maximum route duration and capacity of the trucks available at this depot
For each customer: customer ID, coordinates x and y, service time, demand
For each depot: depot ID, coordinates x and y.

Program

The Hexaly model for the Multi Depot Vehicle Routing Problem (MDVRP) uses list decision variables. For each depot and each truck available at this depot, 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.

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 constrain this quantity to be lower than the truck capacity.

We then compute the distance traveled by each truck. Using another lambda function, we sum up the distances from one customer to the next along the route. Similarly, we compute the total duration of each route and constrain it to be lower than the truck’s maximum route duration.

The objective is to minimize the total distance traveled by all trucks.

Hexaly Modeler
Python
C++
C#
Java

Download

Open in the studio

Execution: hexaly mdvrp.hxm inFileName=instances/p01 [hxTimeLimit=] [outputFileName=]

use io;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

// Compute the distance matrix for customers and for depots/warehouses
function computeDistanceMatrix() {
    for [i in 0...nbCustomers] {
        distanceMatrix[i][i] = 0;
        for [j in (i + 1)...nbCustomers] {
            local localDistance = computeDist(i, j);
            distanceMatrix[j][i] = localDistance;
            distanceMatrix[i][j] = localDistance;
        }
    }

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

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

function computeDistDepot(d, j) {
    return sqrt(pow(depotX[d] - nodesX[j], 2) + pow(depotY[d] - nodesY[j], 2));
}

Execution (Windows): set PYTHONPATH=%HX_HOME%\bin\python

python mdvrp.py instances/p01

Execution (Linux): export PYTHONPATH=/opt/hexaly_13_0/bin/python

python mdvrp.py instances/p01

import hexaly.optimizer
import sys
import math


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

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

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

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

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

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

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

        # Total distance traveled
        total_distance = model.sum()

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

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

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

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

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

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

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

        model.close()

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

        optimizer.solve()

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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


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

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

    return distance_warehouse


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


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

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

    main(instance_file, str_time_limit, output_file)

Compilation / Execution (Windows): cl /EHsc mdvrp.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib

mdvrp instances/p01

Compilation / Execution (Linux): g++ mdvrp.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o mdvrp

./mdvrp instances/p01

#include "optimizer/hexalyoptimizer.h"

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

using namespace hexaly;
using namespace std;

class Mdvrp {
public:
    HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    // Number of trucks per depot
    int nbTrucksPerDepot;

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

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

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

    // Demand per customer
    vector<int> demandsData;

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

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

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

    // Distance traveled by all the trucks
    HxExpression totalDistance;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

        infile.close();
    }

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

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

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

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

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

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

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

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

    return 0;
}

Compilation / Execution (Windows): copy %HX_HOME%\bin\Hexaly.NET.dll .

csc Mdvrp.cs /reference:Hexaly.NET.dll

Mdvrp instances\p01

using Hexaly.Optimizer;
using System;
using System.IO;
using System.Collections.Generic;

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

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    //Number of trucks per depot
    int nbTrucksPerDepot;

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

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

    // Service time per customer
    long[] serviceTimeData;

    // Demand per customer
    long[] demandsData;

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

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

    // Decision variable
    HxExpression[][] customersSequences;

    // Distance traveled by all the trucks
    HxExpression totalDistance;

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

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

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

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

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

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

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

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

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

        HxExpression distMatrix = model.Array(distanceMatrixCustomersData);

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

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

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

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

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

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

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

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

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

        model.Close();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    }

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

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

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

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

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

}

Compilation / Execution (Windows): javac Mdvrp.java -cp %HX_HOME%\bin\hexaly.jar

java -cp %HX_HOME%\bin\hexaly.jar;. Mdvrp instances\p01

Compilation / Execution (Linux): javac Mdvrp.java -cp /opt/hexaly_13_0/bin/hexaly.jar

java -cp /opt/hexaly_13_0/bin/hexaly.jar:. Mdvrp instances/p01

import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

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

    // Number of customers
    int nbCustomers;

    // Number of depots/warehouses
    int nbDepots;

    // Number of trucks per depot
    private int nbTrucksPerDepot;

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

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

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

    // Demand per customer
    private long[] demandsData;

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

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

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

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

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

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

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

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

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

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

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

        HxExpression distMatrix = model.array(distanceMatrixCustomersData);

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

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

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

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

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

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

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

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

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

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

        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}

Problem

Principles learned

Data

Program

Similar examples