Vehicle Routing Problem with Transshipment Facilities (VRPTF)

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_facilities, capacity, customers_demands, \
        depot_distances_data, distance_matrix_data, assignement_costs_data = read_input(instance_file)

    # A point is either a customer or a facility
    # Facilities are duplicated for each customer
    nb_points = nb_customers + nb_customers * nb_facilities

    demands_data = [None] * nb_points
    for c in range(nb_customers):
        demands_data[c] = customers_demands[c]
        for f in range(nb_facilities):
            demands_data[nb_customers + c * nb_facilities + f] = customers_demands[c]

    min_nb_trucks = int(math.ceil(sum(customers_demands) / capacity))
    nb_trucks = int(math.ceil(1.5 * min_nb_trucks))

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

        # Each route is represented as a list containing the points in the order they are visited
        routes_sequences = [m.list(nb_points) for _ in range(nb_trucks)]
        routes = m.array(routes_sequences)

        # Each point must be visited at most once
        m.constraint(m.disjoint(routes_sequences))

        dist_routes = [None] * nb_trucks
        assignement_cost_routes = [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(depot_distances_data)
        assignement_costs = m.array(assignement_costs_data)
        for i in range(nb_points):
            dist_matrix.add_operand(m.array(distance_matrix_data[i]))

        for c in range(nb_customers):
            start_facilities = nb_customers + c * nb_facilities
            end_facilities = start_facilities + nb_facilities

            # Each customer is either contained in a route or assigned to a facility
            facility_used = [m.contains(routes, f) for f in range(start_facilities, end_facilities)]
            delivery_count = m.contains(routes, c) + m.sum(facility_used)
            m.constraint(delivery_count == 1)

        for r in range(nb_trucks):
            route = routes_sequences[r]
            c = m.count(route)

            # Each truck cannot carry more than its capacity
            demand_lambda = m.lambda_function(lambda j: demands[j])
            quantity_served = m.sum(route, demand_lambda)
            m.constraint(quantity_served <= capacity)

            # Distance traveled by each truck
            dist_lambda = m.lambda_function(
                lambda i: m.at(dist_matrix, route[i], route[i + 1]))
            dist_routes[r] = m.sum(m.range(0, c - 1), dist_lambda) + m.iif(
                c > 0,
                m.at(dist_depots, route[0])
                + m.at(dist_depots, route[c - 1]),
                0)
            
            # Cost to assign customers to their facility
            assignment_cost_lambda = m.lambda_function(
                lambda i: assignement_costs[i]
            )
            assignement_cost_routes[r] = m.sum(route, assignment_cost_lambda)

        # The total distance travelled
        total_distance_cost = m.sum(dist_routes)
        # The total assignement cost
        total_assignement_cost = m.sum(assignement_cost_routes)

        # Objective: minimize the sum of the total distance travelled and the total assignement cost
        total_cost = total_distance_cost + total_assignement_cost

        m.minimize(total_cost)

        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: {}; totalCost = {}; totalDistance = {}; totalAssignementCost = {}\n"
                           .format(instance_file, total_cost.value, total_distance_cost.value, total_assignement_cost.value))
                for r in range(nb_trucks):
                    route = routes_sequences[r].value
                    if len(route) == 0:
                        continue
                    file.write("Route {} [".format(r))
                    for i, point in enumerate(route):
                        if point < nb_customers:
                            file.write("Customer {}".format(point))
                        else:
                            file.write("Facility {} assigned to Customer {}"
                                       .format(point % nb_customers, (point - nb_customers) // nb_facilities))
                        if i < len(route) - 1:
                            file.write(", ")
                    file.write("]\n")


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

    nb_customers = int(next(file_it))
    nb_facilities = int(next(file_it))

    facilities_x = [None] * nb_facilities
    facilities_y = [None] * nb_facilities
    for i in range(nb_facilities):
        facilities_x[i] = int(next(file_it))
        facilities_y[i] = int(next(file_it))

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

    truck_capacity = int(next(file_it))

    # Facility capacities : skip
    for f in range(nb_facilities):
        next(file_it)

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

    depot_x, depot_y = compute_depot_coordinates(customers_x, customers_y,
                                                 facilities_x, facilities_y)
    depot_distances, distance_matrix = compute_distances(customers_x, customers_y,
                                                         facilities_x, facilities_y,
                                                         depot_x, depot_y)
    assignement_costs = compute_assignment_costs(nb_customers, nb_facilities, distance_matrix)

    return nb_customers, nb_facilities, truck_capacity, customer_demands, \
        depot_distances, distance_matrix, assignement_costs


def compute_depot_coordinates(customers_x, customers_y, facilities_x, facilities_y):
    # Compute the coordinates of the bounding box containing all of the points
    x_min = min(min(customers_x), min(facilities_x))
    x_max = max(max(customers_x), max(facilities_x))
    y_min = min(min(customers_y), min(facilities_y))
    y_max = max(max(customers_y), max(facilities_y))

    # We assume that the depot is at the center of the bounding box
    return x_min + (x_max - x_min) // 2, y_min + (y_max - y_min) // 2


def compute_distances(customers_x, customers_y, facilities_x, facilities_y, depot_x, depot_y):
    nb_customers = len(customers_x)
    nb_facilities = len(facilities_x)
    nb_points = nb_customers + nb_customers * nb_facilities


    # Distance to depot
    depot_distances = [None] * nb_points

    # Customer to depot
    for c in range(nb_customers):
        depot_distances[c] = compute_dist(customers_x[c], depot_x, customers_y[c], depot_y)

    # Facility to depot
    for c in range(nb_customers):
        for f in range(nb_facilities):
            depot_distances[nb_customers + c * nb_facilities + f] = \
                compute_dist(facilities_x[f], depot_x, facilities_y[f], depot_y)

    # Distance between points
    distance_matrix = [[None for _ in range(nb_points)] for _ in range(nb_points)]

    # Distances between customers
    for c_1 in range(nb_customers):
        for c_2 in range(nb_customers):
            distance_matrix[c_1][c_2] = \
                compute_dist(customers_x[c_1], customers_x[c_2],
                             customers_y[c_1], customers_y[c_2])

    # Distances between customers and facilities
    for c_1 in range(nb_customers):
        for f in range(nb_facilities):
            distance = compute_dist(facilities_x[f], customers_x[c_1],
                                    facilities_y[f], customers_y[c_1])
            for c_2 in range(nb_customers):
                # Index representing serving c_2 through facility f
                facility_index = nb_customers + c_2 * nb_facilities + f
                distance_matrix[facility_index][c_1] = distance
                distance_matrix[c_1][facility_index] = distance

    # Distances between facilities
    for f_1 in range(nb_facilities):
        for f_2 in range(nb_facilities):
            dist = compute_dist(facilities_x[f_1], facilities_x[f_2], facilities_y[f_1], facilities_y[f_2])
            for c_1 in range(nb_customers):
                for c_2 in range(nb_customers):
                    distance_matrix[nb_customers + c_1 * nb_facilities + f_1]\
                        [nb_customers + c_2 * nb_facilities + f_2] = dist

    return depot_distances, distance_matrix


def compute_assignment_costs(nb_customers, nb_facilities, distance_matrix):
    # Compute assignment cost for each point
    nb_points = nb_customers + nb_customers * nb_facilities
    assignment_costs = [0] * nb_points
    for c in range(nb_customers):
        for f in range(nb_facilities):
            #  Cost of serving customer c through facility f
            assignment_costs[nb_customers + c * nb_facilities + f] = \
                distance_matrix[c][nb_customers + c * nb_facilities + f]
    return assignment_costs


def compute_dist(xi, xj, yi, yj):
    exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
    return round(exact_dist)


def read_input(filename):
    if filename.endswith(".dat"):
        return read_input_dat(filename)
    else:
        raise Exception("Unknown file format")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python vrptf.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)

#include <climits>
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <numeric>
#include <vector>
#include <algorithm>

#include "optimizer/hexalyoptimizer.h"

using namespace hexaly;
using namespace std;

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

    // Number of customers
    int nbCustomers;

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

    // Customers demands
    vector<int> demandsData;

    // Number of depots
    int nbFacilities;

    // Depots coordinates
    vector<int> xFacilities;
    vector<int> yFacilities;

    // Number of points
    int nbPoints;

    // Depot coordinates
    int xDepot;
    int yDepot;

    // Number of trucks
    int nbTrucks;

    // Capacity of trucks
    int truckCapacity;

    // Distance matrix
    vector<vector<int>> distMatrixData;

    // Distance to depot array
    vector<int> distDepotsData;

    // Assignement costs
    vector<int> assignmentCostsData;

    // Decision variables
    vector<HxExpression> routesSequences;

    // Objective value
    HxExpression totalDistanceCost;
    HxExpression totalAssignementCost;
    HxExpression totalCost;

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

    void solve(const char* limit) {
        // Declare the optimization model
        HxModel m = optimizer.getModel();

        double totalDemand = accumulate(demandsData.begin(), demandsData.begin() + nbCustomers, 0);
        int minNbTrucks = ceil(totalDemand / truckCapacity);
        nbTrucks = ceil(1.5 * minNbTrucks);

        // A route is represented as a list containing the points in the order they are visited
        routesSequences.resize(nbTrucks);
        for (int r = 0; r < nbTrucks; ++r) {
            routesSequences[r] = m.listVar(nbPoints);
        }

        HxExpression routes = m.array(routesSequences.begin(), routesSequences.end());

        // Each customer must be visited at most once
        m.constraint(m.disjoint(routesSequences.begin(), routesSequences.end()));

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

        vector<HxExpression> distRoutes(nbTrucks);
        vector<HxExpression> assignmentCostRoutes(nbTrucks);

        for (int c = 0; c < nbCustomers; ++c) {
            int startFacilities = nbCustomers + c * nbFacilities;
            int endFacilities = startFacilities + nbFacilities;

            // Each customer is either contained in a route or assigned to a facility
            HxExpression facilityUsedSum = m.sum();
            for (int f = startFacilities; f < endFacilities; ++f) {
                facilityUsedSum.addOperand(m.contains(routes, f));
            }
            m.constraint(m.contains(routes, c) + facilityUsedSum == 1);
        }

        for (int r = 0; r < nbTrucks; ++r) {
            HxExpression route = routesSequences[r];
            HxExpression c = m.count(route);

            // Each truck cannot carry more than its capacity
            HxExpression demandLambda = m.lambdaFunction([&](HxExpression j) { return demands[j]; });
            HxExpression quantityServed = m.sum(route, demandLambda);
            m.constraint(quantityServed <= truckCapacity);

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

            // Truck is used if it visits at least one point
            HxExpression truckUsed = c > 0;

            // Distance traveled by each truck
            distRoutes[r] = m.sum(m.range(0, c - 1), distLambda) + 
                m.iif(truckUsed,
                    m.at(distDepots, route[0]) +
                    m.at(distDepots, route[c - 1]),
                    0);

            // The cost to assign customers to their facility
            HxExpression assignementCostLambda =
                m.lambdaFunction([&](HxExpression i) { return assignementCosts[i]; });
            assignmentCostRoutes[r] = m.sum(route, assignementCostLambda);
        }

        // The total distance travelled
        totalDistanceCost = m.sum(distRoutes.begin(), distRoutes.end());
        // The total assignement cost
        totalAssignementCost = m.sum(assignmentCostRoutes.begin(), assignmentCostRoutes.end());

        // Objective: minimize the sum of the total distance travelled and the total assignement cost
        totalCost = totalDistanceCost + totalAssignementCost;

        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 << "; totalCost = " << totalCost.getIntValue()
            << "; totalDistance = " << totalDistanceCost.getIntValue()
            << "; totalAssignementCost = " << totalAssignementCost.getIntValue() << endl;
        for (int r = 0; r < nbTrucks; ++r) {
            HxCollection route = routesSequences[r].getCollectionValue();
            if (route.count() == 0) continue;
            file << "Route " << r << " [";
            for (int i = 0; i < route.count(); ++i) {
                long point = route.get(i);
                if (point < nbCustomers) {
                    file << "Customer " << point;
                }
                else {
                    file << "Facility " << point % nbCustomers << " assigned to Customer " << (point - nbCustomers) / nbFacilities;
                }
                if (i < route.count() - 1) {
                    file << ", ";
                }
            }
            file << "]" << endl;
        }
    }

private:
    void readInput(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);
        infile >> nbFacilities;
        xFacilities.resize(nbFacilities);
        yFacilities.resize(nbFacilities);

        // A point is either a customer or a facility
        // Facilities are duplicated for each customer
        nbPoints = nbCustomers + nbCustomers * nbFacilities;
        demandsData.resize(nbPoints);
        distMatrixData.resize(nbPoints);
        for (int i = 0; i < nbPoints; ++i) {
            distMatrixData[i].resize(nbPoints);
        }
        distDepotsData.resize(nbPoints);

        for (int f = 0; f < nbFacilities; ++f) {
            infile >> xFacilities[f];
            infile >> yFacilities[f];
        }
        for (int c = 0; c < nbCustomers; ++c) {
            infile >> xCustomers[c];
            infile >> yCustomers[c];
        }

        infile >> truckCapacity;

        // Facility capacities : skip
        infile.ignore(1000, '\n');
        infile.ignore(1000, '\n');
        for (int f = 0; f < nbFacilities; ++f) {
            infile.ignore(1000, '\n');
        }

        for (int c = 0; c < nbCustomers; ++c) {
            infile >> demandsData[c];
            for (int f = 0; f < nbFacilities; ++f) {
                demandsData[nbCustomers + c * nbFacilities + f] = demandsData[c];
            }
        }

        computeDepotCoordinates();
        computeDistances();
        computeAssignmentCosts();
    }

    void computeDepotCoordinates() {
        // Compute the coordinates of the bounding box containing all of the points
        int xMin = INT_MAX;
        int xMax = INT_MIN;
        int yMin = INT_MAX;
        int yMax = INT_MIN;

        for (int c = 0; c < nbCustomers; ++c) {
            xMin = min(xMin, xCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            xMin = min(xMin, xFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            xMax = max(xMax, xCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            xMax = max(xMax, xFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            yMin = min(yMin, yCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            yMin = min(yMin, yFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            yMax = max(yMax, yCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            yMax = max(yMax, yFacilities[f]);
        }

        // We assume that the depot is at the center of the bounding box
        xDepot = xMin + (xMax - xMin) / 2;
        yDepot = yMin + (yMax - yMin) / 2;
    }

    void computeDistances() {
        // Customer to depot
        for (int c = 0; c < nbCustomers; ++c) {
            distDepotsData[c] = computeDist(xCustomers[c], xDepot, yCustomers[c], yDepot);
        }

        // Facility to depot
        for (int c = 0; c < nbCustomers; ++c) {
            for (int f = 0; f < nbFacilities; ++f) {
                distDepotsData[nbCustomers + c * nbFacilities + f] = computeDist(xFacilities[f], xDepot, yFacilities[f],
                    yDepot);
            }
        }

        // Distances between customers
        for (int c1 = 0; c1 < nbCustomers; ++c1) {
            for (int c2 = 0; c2 < nbCustomers; ++c2) {
                long dist = computeDist(xCustomers[c1], xCustomers[c2], yCustomers[c1], yCustomers[c2]);
                distMatrixData[c1][c2] = dist;
                distMatrixData[c2][c1] = dist;
            }
        }

        // Distances between customers and facilities
        for (int c1 = 0; c1 < nbCustomers; ++c1) {
            for (int f = 0; f < nbFacilities; ++f) {
                long dist = computeDist(xFacilities[f], xCustomers[c1], yFacilities[f], yCustomers[c1]);
                for (int c2 = 0; c2 < nbCustomers; ++c2) {
                     // Index representing serving c2 through facility f
                    int facilityIndex = nbCustomers + c2 * nbFacilities + f;
                    distMatrixData[facilityIndex][c1] = dist;
                    distMatrixData[c1][facilityIndex] = dist;
                }
            }
        }

        // Distances between facilities
        for (int f1 = 0; f1 < nbFacilities; ++f1) {
            for (int f2 = 0; f2 < nbFacilities; ++f2) {
                long dist = computeDist(xFacilities[f1], xFacilities[f2], yFacilities[f1], yFacilities[f2]);
                for (int c1 = 0; c1 < nbCustomers; ++c1) {
                    // Index representing serving c1 through facility f1
                    int index1 = nbCustomers + c1 * nbFacilities + f1;
                    for (int c2 = 0; c2 < nbCustomers; ++c2) {
                        // Index representing serving c2 through facility f2
                        int index2 = nbCustomers + c2 * nbFacilities + f2;
                        distMatrixData[index1][index2] = dist;
                    }
                }
            }
        }
    }

    void computeAssignmentCosts() {
        // Compute assignment cost for each point
        assignmentCostsData.resize(nbPoints);
        for (int c = 0; c < nbCustomers; ++c) {
            assignmentCostsData[c] = 0;
            for (int f = 0; f < nbFacilities; ++f) {
                // Cost of serving customer c through facility f
                assignmentCostsData[nbCustomers + c * nbFacilities + f] =
                    distMatrixData[c][nbCustomers + c * nbFacilities + f];
            }
        }
    }

    int computeDist(int xi, int xj, int yi, int yj) {
        double exactDist = sqrt(pow(xi - xj, 2) + pow(yi - yj, 2));
        return round(exactDist);
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: ./vrptf 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 {
        Vrptf 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;
}

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

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

    // Number of customers
    int nbCustomers;

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

    // Customers demands
    int[] demandsData;

    // Number of facilities
    int nbFacilities;

    // Depots coordinates
    int[] xFacilities;
    int[] yFacilities;

    // Number of points
    int nbPoints;

    // Depot coordinates
    int xDepot;
    int yDepot;

    // Number of trucks
    int nbTrucks;

    // Capacity of a truck
    int truckCapacity;

    // Distance matrix
    long[][] distMatrixData;

    // Distance to depot array
    long[] distDepotsData;

    // Assignement costs
    long[] assignmentCostsData;

    // Decision variables
    HxExpression[] routesSequences;

    // Objective value
    HxExpression totalDistanceCost;
    HxExpression totalAssignementCost;
    HxExpression totalCost;

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

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

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

        double totalDemand = 0;
        for (int c = 0; c < nbCustomers; ++c)
        {
            totalDemand += demandsData[c];
        }
        int minNbTrucks = (int)Math.Ceiling(totalDemand / truckCapacity);
        nbTrucks = (int)Math.Ceiling(1.5 * minNbTrucks);

        routesSequences = new HxExpression[nbTrucks];

        // A route is represented as a list containing the points in the order they are visited
        for (int r = 0; r < nbTrucks; ++r)
        {
            routesSequences[r] = m.List(nbPoints);
        }

        HxExpression routes = m.Array(routesSequences);

        // Each point must be visited at most once
        m.Constraint(m.Disjoint(routesSequences));

        HxExpression demands = m.Array(demandsData);
        HxExpression distMatrix = m.Array(distMatrixData);
        HxExpression distDepots = m.Array(distDepotsData);
        HxExpression assignmentCosts = m.Array(assignmentCostsData);

        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] assignmentCostRoutes = new HxExpression[nbTrucks];

        for (int c = 0; c < nbCustomers; ++c)
        {
            int startFacilities = nbCustomers + c * nbFacilities;
            int endFacilities = startFacilities + nbFacilities;

            // Each customer is either contained in a route or assigned to a facility
            HxExpression facilityUsedSum = m.Sum();
            for (int f = startFacilities; f < endFacilities; ++f)
            {
                facilityUsedSum.AddOperand(m.Contains(routes, f));
            }
            m.Constraint(m.Contains(routes, c) + facilityUsedSum == 1);
        }

        for (int r = 0; r < nbTrucks; ++r)
        {
            HxExpression route = routesSequences[r];
            HxExpression c = m.Count(route);

            // Each truck cannot carry more than its capacity
            HxExpression demandLambda = m.LambdaFunction(j => m.At(demands, j));
            HxExpression quantityServed = m.Sum(route, demandLambda);
            m.Constraint(quantityServed <= truckCapacity);

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

            // Truck is used if it visits at least one point
            HxExpression truckUsed = c > 0;

            // Distance traveled by each truck
            distRoutes[r] = m.Sum(m.Range(0, c - 1), distLambda) +
                    m.If(truckUsed,
                    m.At(distDepots, m.At(route, 0)) +
                            m.At(distDepots, m.At(route, c - 1)),
                            0);

            // The cost to assign customers to their facility
            HxExpression assignmentCostLambda = m.LambdaFunction(i => m.At(assignmentCosts, i));
            assignmentCostRoutes[r] = m.Sum(route, assignmentCostLambda);
        }

        // The total distance travelled
        totalDistanceCost = m.Sum(distRoutes);
        // The total assignement cost
        totalAssignementCost = m.Sum(assignmentCostRoutes);

        // Objective: minimize the sum of the total distance travelled and the total assignement cost
        totalCost = totalDistanceCost + totalAssignementCost;

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

        optimizer.GetParam().SetTimeLimit(limit);
        optimizer.Solve();
    }

    /* 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];

            line = ReadNextLine(input);
            nbFacilities = int.Parse(line[0]);
            xFacilities = new int[nbFacilities];
            yFacilities = new int[nbFacilities];

            // A point is either a customer or a facility
            // Facilities are duplicated for each customer
            nbPoints = nbCustomers + nbCustomers * nbFacilities;
            distMatrixData = new long[nbPoints][];
            distDepotsData = new long[nbPoints];
            demandsData = new int[nbPoints];

            line = ReadNextLine(input);
            for (int f = 0; f < nbFacilities; ++f)
            {
                line = ReadNextLine(input);
                xFacilities[f] = int.Parse(line[0]);
                yFacilities[f] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            for (int c = 0; c < nbCustomers; ++c)
            {
                line = ReadNextLine(input);
                xCustomers[c] = int.Parse(line[0]);
                yCustomers[c] = int.Parse(line[1]);
            }
            line = ReadNextLine(input);
            line = ReadNextLine(input);
            truckCapacity = int.Parse(line[0]);
            line = ReadNextLine(input);

            // Facility capacities : skip
            for (int f = 0; f < nbFacilities; ++f)
            {
                ReadNextLine(input);
            }
            line = ReadNextLine(input);
            for (int c = 0; c < nbCustomers; ++c)
            {
                line = ReadNextLine(input);
                demandsData[c] = int.Parse(line[0]);
                for (int f = 0; f < nbFacilities; f++)
                {
                    demandsData[nbCustomers + c * nbFacilities + f] = demandsData[c];
                }
            }
            line = ReadNextLine(input);

            ComputeDepotCoordinates();
            ComputeDistanceArrays();
            ComputeAssignmentCosts();
        }
    }

    void ComputeDepotCoordinates()
    {
        // Compute the coordinates of the bounding box containing all of the points
        int xMin = Math.Min(xCustomers.Min(), xFacilities.Min());
        int xMax = Math.Max(xCustomers.Max(), xFacilities.Max());
        int yMin = Math.Min(yCustomers.Min(), yFacilities.Min());
        int yMax = Math.Max(yCustomers.Max(), yFacilities.Max());

        // We assume that the depot is at the center of the bounding box
        xDepot = xMin + (xMax - xMin) / 2;
        yDepot = yMin + (yMax - yMin) / 2;
    }

    void ComputeDistanceArrays()
    {
        distDepotsData = new long[nbPoints];

        // Customer to depot
        for (int c = 0; c < nbCustomers; ++c)
        {
            distDepotsData[c] = ComputeDist(xCustomers[c], xDepot, yCustomers[c], yDepot);
        }

        // Facility to depot
        for (int c = 0; c < nbCustomers; ++c)
        {
            for (int f = 0; f < nbFacilities; ++f)
            {
                distDepotsData[nbCustomers + c * nbFacilities + f] = ComputeDist(xFacilities[f], xDepot, yFacilities[f], yDepot);
            }
        }

        distMatrixData = new long[nbPoints][];
        for (int i = 0; i < nbPoints; i++)
        {
            distMatrixData[i] = new long[nbPoints];
        }

        // Distances between customers
        for (int c1 = 0; c1 < nbCustomers; ++c1)
        {
            for (int c2 = 0; c2 < nbCustomers; ++c2)
            {
                long dist = ComputeDist(xCustomers[c1], xCustomers[c2], yCustomers[c1], yCustomers[c2]);
                distMatrixData[c1][c2] = dist;
                distMatrixData[c2][c1] = dist;
            }
        }

        // Distances between customers and facilities
        for (int c1 = 0; c1 < nbCustomers; ++c1)
        {
            for (int f = 0; f < nbFacilities; ++f)
            {
                long dist = ComputeDist(xFacilities[f], xCustomers[c1], yFacilities[f], yCustomers[c1]);
                for (int c2 = 0; c2 < nbCustomers; ++c2)
                {
                    // Index representing serving c2 through facility f
                    int facilityIndex = nbCustomers + c2 * nbFacilities + f;
                    distMatrixData[facilityIndex][c1] = dist;
                    distMatrixData[c1][facilityIndex] = dist;
                }
            }
        }
    
        // Distances between facilities
        for (int f1 = 0; f1 < nbFacilities; ++f1)
        {
            for (int f2 = 0; f2 < nbFacilities; ++f2)
            {
                long dist = ComputeDist(xFacilities[f1], xFacilities[f2], yFacilities[f1], yFacilities[f2]);
                for (int c1 = 0; c1 < nbCustomers; ++c1)
                {
                    // Index representing serving c1 through facility f1
                    int index1 = nbCustomers + c1 * nbFacilities + f1;
                    for (int c2 = 0; c2 < nbCustomers; ++c2)
                    {
                        // Index representing serving c2 through facility f2
                        int index2 = nbCustomers + c2 * nbFacilities + f2;
                        distMatrixData[index1][index2] = dist;
                    }
                }
            }
        }
    }

    private long ComputeDist(int xi, int xj, int yi, int yj)
    {
        double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        return Convert.ToInt64(Math.Round(exactDist));
    }

    private void ComputeAssignmentCosts() {
        // Compute assignment cost for each point
        assignmentCostsData = new long[nbPoints];
        for (int c = 0; c < nbCustomers; ++c)
        {
            assignmentCostsData[c] = 0;
            for (int f = 0; f < nbFacilities; ++f)
            {
                // Cost of serving customer c through facility f
                assignmentCostsData[nbCustomers + c * nbFacilities + f] =
                    distMatrixData[c][nbCustomers + c * nbFacilities + f];
            }
        }
    }

    /* Write the solution in a file */
    void WriteSolution(string infile, string outfile)
    {
        using (StreamWriter output = new StreamWriter(outfile))
        {
            output.WriteLine(
                "File name: " + infile + "; totalCost = " + totalCost.GetIntValue()
                + "; totalDistance = " + totalDistanceCost.GetIntValue()
                + "; totalAssignementCost = " + totalAssignementCost.GetIntValue()
            );
            for (int r = 0; r < nbTrucks; ++r)
            {
                HxCollection route = routesSequences[r].GetCollectionValue();
                if (route.Count() == 0) continue;
                output.Write("Route " + r + " [");
                for (int i = 0; i < route.Count(); i++)
                {
                    long point = route.Get(i);
                    if (point < nbCustomers)
                    {
                        output.Write("Customer " + point);
                    }
                    else
                    {
                        output.Write("Facility " + point % nbCustomers + " assigned to Customer " + (point - nbCustomers) / nbFacilities);
                    }
                    if (i < route.Count() - 1)
                    {
                        output.Write(", ");
                    }
                }
                output.WriteLine("]");
            }
        }
    }

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

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Vrptf 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 (Vrptf model = new Vrptf())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (strOutput != null)
                model.WriteSolution(instanceFile, strOutput);
        }
    }
}

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

public class Vrptf {
    // 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 facilities
    private int nbFacilities;

    // Facilities coordinates
    private int[] xFacilities;
    private int[] yFacilities;

    // Number of points
    private int nbPoints;

    // Depot coordinates
    private int xDepot;
    private int yDepot;

    // Number of trucks
    private int nbTrucks;

    // Capacity of trucks
    private int truckCapacity;

    // Distance matrix
    private long[][] distMatrixData;

    // Distance to depot array
    private long[] distDepotsData;

    // Assignement costs
    private long[] assignmentCostsData;

    // Decision variables
    private HxExpression[] routesSequences;

    // Objective value
    private HxExpression totalDistanceCost;
    private HxExpression totalAssignementCost;
    private HxExpression totalCost;

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

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

        double totalDemand = 0;
        for (int c = 0; c < nbCustomers; ++c) {
            totalDemand += demandsData[c];
        }
        int minNbTrucks = (int) Math.ceil(totalDemand / truckCapacity);
        nbTrucks = (int) Math.ceil(1.5 * minNbTrucks);

        routesSequences = new HxExpression[nbTrucks];

        // A route is represented as a list containing the points in the order they are
        // visited
        for (int i = 0; i < nbTrucks; ++i) {
            routesSequences[i] = m.listVar(nbPoints);
        }

        HxExpression routes = m.array(routesSequences);

        // Each point must be visited at most once
        m.constraint(m.disjoint(routesSequences));

        HxExpression demands = m.array(demandsData);
        HxExpression distMatrix = m.array(distMatrixData);
        HxExpression distDepots = m.array(distDepotsData);
        HxExpression assignmentCosts = m.array(assignmentCostsData);

        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] assignmentCostRoutes = new HxExpression[nbTrucks];

        for (int c = 0; c < nbCustomers; ++c) {
            int startFacilities = nbCustomers + c * nbFacilities;
            int endFacilities = startFacilities + nbFacilities;

            // Each customer is either contained in a route or assigned to a facility
            HxExpression facilityUsedSum = m.sum();
            for (int f = startFacilities; f < endFacilities; ++f) {
                facilityUsedSum.addOperand(m.contains(routes, f));
            }
            m.constraint(m.eq(m.sum(m.contains(routes, c), facilityUsedSum), 1));
        }

        for (int r = 0; r < nbTrucks; ++r) {
            HxExpression route = routesSequences[r];
            HxExpression c = m.count(route);

            // Each truck cannot carry more than its capacity
            HxExpression demandLambda = m.lambdaFunction(j -> m.at(demands, j));
            HxExpression quantityServed = m.sum(route, demandLambda);
            m.constraint(m.leq(quantityServed, truckCapacity));

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

            // Truck is used if it visits at least one point
            HxExpression truckUsed = m.gt(c, 0);

            // Distance traveled by each truck
            distRoutes[r] = m.sum(
                    m.sum(m.range(0, m.sub(c, 1)), distLambda),
                    m.iif(truckUsed,
                            m.sum(m.at(distDepots, m.at(route, 0)),
                                    m.at(distDepots, m.at(route, m.sub(c, 1)))),
                            0));

            // The cost to assign customers to their facility
            HxExpression assignmentCostLambda = m.lambdaFunction(i -> m.at(assignmentCosts, i));
            assignmentCostRoutes[r] = m.sum(route, assignmentCostLambda);
        }

        // The total distance travelled
        totalDistanceCost = m.sum(distRoutes);
        // The total assignement cost
        totalAssignementCost = m.sum(assignmentCostRoutes);

        // Objective: minimize the sum of the total distance travelled and the total
        // assignement cost
        totalCost = m.sum(totalDistanceCost, totalAssignementCost);
        m.minimize(totalCost);
        m.close();

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

    /* Write the solution to a file */
    private void writeSolution(String infile, String outfile) throws IOException {
        try (PrintWriter output = new PrintWriter(outfile)) {
            output.println("File name: " + infile + "; totalCost = " + totalCost.getIntValue() + "; totalDistance = "
                    + totalDistanceCost.getIntValue() + "; totalAssignementCost = "
                    + totalAssignementCost.getIntValue());
            for (int r = 0; r < nbTrucks; ++r) {
                HxCollection route = routesSequences[r].getCollectionValue();
                if (route.count() == 0)
                    continue;
                output.print("Route " + r + " [");
                for (int i = 0; i < route.count(); i++) {
                    long point = route.get(i);
                    if (point < nbCustomers) {
                        output.print("Customer " + point);
                    } else {
                        output.print("Facility " + point % nbCustomers + " assigned to Customer "
                                + (point - nbCustomers) / nbFacilities);
                    }
                    if (i < route.count() - 1) {
                        output.print(", ");
                    }
                }
                output.println("]");
            }
        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.useLocale(Locale.ROOT);
            nbCustomers = input.nextInt();
            nbFacilities = input.nextInt();
            // A point is either a customer or a facility
            // Facilities are duplicated for each customer
            nbPoints = nbCustomers + nbCustomers * nbFacilities;

            xFacilities = new int[nbFacilities];
            yFacilities = new int[nbFacilities];
            for (int f = 0; f < nbFacilities; ++f) {
                xFacilities[f] = input.nextInt();
                yFacilities[f] = input.nextInt();
            }

            xCustomers = new int[nbCustomers];
            yCustomers = new int[nbCustomers];
            for (int c = 0; c < nbCustomers; ++c) {
                xCustomers[c] = input.nextInt();
                yCustomers[c] = input.nextInt();
            }

            truckCapacity = input.nextInt();

            // Facility capacities : skip
            for (int f = 0; f < nbFacilities; ++f) {
                input.next();
            }

            demandsData = new int[nbPoints];
            for (int c = 0; c < nbCustomers; ++c) {
                demandsData[c] = input.nextInt();
                for (int f = 0; f < nbFacilities; ++f) {
                    demandsData[nbCustomers + c * nbFacilities + f] = demandsData[c];
                }
            }

            computeDepotCoordinates();
            computeDistances();
            computeAssignmentCosts();
        } catch (FileNotFoundException e) {
            System.out.println("An error occurred.");
            e.printStackTrace();
        }
    }

    void computeDepotCoordinates() {
        // Compute the coordinates of the bounding box containing all of the points
        int xMin = Integer.MAX_VALUE;
        int xMax = Integer.MIN_VALUE;
        int yMin = Integer.MAX_VALUE;
        int yMax = Integer.MIN_VALUE;

        for (int c = 0; c < nbCustomers; ++c) {
            xMin = Math.min(xMin, xCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            xMin = Math.min(xMin, xFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            xMax = Math.max(xMax, xCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            xMax = Math.max(xMax, xFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            yMin = Math.min(yMin, yCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            yMin = Math.min(yMin, yFacilities[f]);
        }
        for (int c = 0; c < nbCustomers; ++c) {
            yMax = Math.max(yMax, yCustomers[c]);
        }
        for (int f = 0; f < nbFacilities; ++f) {
            yMax = Math.max(yMax, yFacilities[f]);
        }

        // We assume that the depot is at the center of the bounding box
        xDepot = xMin + (xMax - xMin) / 2;
        yDepot = yMin + (yMax - yMin) / 2;
    }

    void computeDistances() {
        distDepotsData = new long[nbPoints];

        // Customer to depot
        for (int c = 0; c < nbCustomers; ++c) {
            distDepotsData[c] = computeDist(xCustomers[c], xDepot, yCustomers[c], yDepot);
        }

        // Facility to depot
        for (int c = 0; c < nbCustomers; ++c) {
            for (int f = 0; f < nbFacilities; ++f) {
                distDepotsData[nbCustomers + c * nbFacilities + f] = computeDist(xFacilities[f], xDepot, yFacilities[f],
                        yDepot);
            }
        }

        distMatrixData = new long[nbPoints][];
        for (int i = 0; i < nbPoints; i++) {
            distMatrixData[i] = new long[nbPoints];
        }

        // Distances between customers
        for (int c1 = 0; c1 < nbCustomers; ++c1) {
            for (int c2 = 0; c2 < nbCustomers; ++c2) {
                long dist = computeDist(xCustomers[c1], xCustomers[c2], yCustomers[c1], yCustomers[c2]);
                distMatrixData[c1][c2] = dist;
                distMatrixData[c2][c1] = dist;
            }
        }

        // Distances between customers and facilities
        for (int c1 = 0; c1 < nbCustomers; ++c1) {
            for (int f = 0; f < nbFacilities; ++f) {
                long dist = computeDist(xFacilities[f], xCustomers[c1], yFacilities[f], yCustomers[c1]);
                for (int c2 = 0; c2 < nbCustomers; ++c2) {
                     // Index representing serving c2 through facility f
                    int facilityIndex = nbCustomers + c2 * nbFacilities + f;
                    distMatrixData[facilityIndex][c1] = dist;
                    distMatrixData[c1][facilityIndex] = dist;
                }
            }
        }

        // Distances between facilities
        for (int f1 = 0; f1 < nbFacilities; ++f1) {
            for (int f2 = 0; f2 < nbFacilities; ++f2) {
                long dist = computeDist(xFacilities[f1], xFacilities[f2], yFacilities[f1], yFacilities[f2]);
                for (int c1 = 0; c1 < nbCustomers; ++c1) {
                    // Index representing serving c1 through facility f1
                    int index1 = nbCustomers + c1 * nbFacilities + f1;
                    for (int c2 = 0; c2 < nbCustomers; ++c2) {
                        // Index representing serving c2 through facility f2
                        int index2 = nbCustomers + c2 * nbFacilities + f2;
                        distMatrixData[index1][index2] = dist;
                    }
                }
            }
        }
    }

    void computeAssignmentCosts() {
        // Compute assignment cost for each point
        assignmentCostsData = new long[nbPoints];
        for (int c = 0; c < nbCustomers; ++c) {
            assignmentCostsData[c] = 0;
            for (int f = 0; f < nbFacilities; ++f) {
                // Cost of serving customer c through facility f
                assignmentCostsData[nbCustomers + c * nbFacilities + f] =
                    distMatrixData[c][nbCustomers + c * nbFacilities + f];
            }
        }
    }

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

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java vrptf 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";

            Vrptf model = new Vrptf(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (strOutfile != null)
                model.writeSolution(instanceFile, strOutfile);
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            throw new RuntimeException(ex);
        }
    }
}