View all examples

Inventory Routing Problem (IRP)

Problem

The Inventory Routing Problem (IRP) is a distribution problem in which a product has to be shipped from a supplier to several customers over a given time horizon. Each customer defines a maximum inventory level. The supplier monitors each customer’s inventory and determines its replenishment policy, guaranteeing that no stockout occurs at the customer (vendor-managed inventory policy). Shipments from the supplier to the customers are performed by a vehicle of a given capacity.

Principles learned

  • Add list decision variables to model the trucks’ sequences of customers
  • Define lambda functions to compute the traveled distance at each time slot
  • Define expressions recursively to compute the inventory levels

Data

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

  • First line: the number of customers, the number of discrete time slots in the planning time horizon, and the vehicle’s transportation capacity.
  • Second line: for the supplier
    • The index of the supplier
    • The x coordinate
    • The y coordinate
    • The starting level of the supplier’s inventory
    • The quantity of product made available at each time slot
    • The unit inventory cost
  • Next lines: for each customer
    • The index of the customer
    • The x coordinate
    • The y coordinate
    • The starting level of the inventory
    • The maximum level of inventory
    • The minimum level of inventory
    • The quantity absorbed at each time slot
    • The unit inventory cost

Program

The Hexaly model for the Inventory Routing Problem (IRP) uses list decision variables. For each discrete time slot, we define a list variable representing the customers visited by the truck during this time slot. The cost of each route depends on the distance traveled by the truck. Using a lambda function, we sum up the distances from one customer to the next along the route.

We compute the supplier’s inventory levels recursively. Indeed, its level at time t is given by its level at time t-1, plus the product quantity made available at time t-1, minus the total quantity shipped to the customers at time t-1. The stockout constraints at the supplier ensure that the supplier’s inventory level is always sufficient to ship the total quantity delivered to the customers during a given time slot.

Similarly, we then compute the customers’ inventory levels. For each customer, the inventory level at time t is given by the level at time t-1, plus the product quantity shipped from the supplier to the customer at time t-1, minus the product quantity consumed at time t-1. The stockout constraints at the customers guarantee that each customer’s inventory level is always positive.

The capacity constraints ensure that the total product quantity loaded on the vehicle never exceeds its capacity. The maximum level constraints guarantee that the product quantity shipped from the supplier to each customer is not greater than the maximum inventory level of each customer.

The objective is to minimize the sum of all costs: the total inventory cost at the supplier, the total inventory cost at customers, and the total transportation cost.

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

/* Read instance data. The input files follow the "Archetti" format */
function input() {
    usage = "Usage: hexaly irp.hxm " + 
            "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;

    readInputIrp();

    // Compute distance matrix
    computeDistanceMatrix();
}

/* Declare the optimization model */
function model() {
    // Quantity of product delivered at each discrete time instant of
    // the planning time horizon to each customer
    delivery[t in 0...horizonLength][i in 0...nbCustomers] <- float(0, capacity);

    // Sequence of customers visited at each discrete time instant of
    // the planning time horizon
    route[t in 0...horizonLength] <- list(nbCustomers);

    for [t in 0...horizonLength] {
        local sequence <- route[t];
        local c <- count(sequence);

        // Customers receive products only if they are visited
        isDelivered[t][i in 0...nbCustomers] <- contains(sequence, i);

        // Distance traveled at instant t
        costRoute[t] <- (c > 0) ? distanceSupplier[sequence[0]] + 
                sum(1...c, i => distance[sequence[i-1]][sequence[i]]) + 
                distanceSupplier[sequence[c-1]] : 0;
    }

    // Stockout constraints at the supplier
    inventorySupplier[0] <- startLevelSupplier;
    for [t in 1...horizonLength+1] {
        inventorySupplier[t] <- inventorySupplier[t-1] 
                - sum[i in 0...nbCustomers](delivery[t-1][i])
                + productionRate_supplier;

        if (t != horizonLength)
            constraint inventorySupplier[t] >= sum[i in 0...nbCustomers](delivery[t][i]);
    }

    // Stockout constraints at the customers
    inventory[i in 0...nbCustomers][0] <- startLevel[i];
    for [i in 0...nbCustomers] {
        for [t in 1...horizonLength+1] {
            inventory[i][t] <- inventory[i][t-1] + delivery[t-1][i] - demandRate[i];
            constraint inventory[i][t] >= 0;
        }
    }
    
    for [t in 0...horizonLength] {
        // Capacity constraints
        constraint sum[i in 0...nbCustomers](delivery[t][i]) <= capacity;
        
        // Maximum level constraints
        for [i in 0...nbCustomers] {
            constraint delivery[t][i] <= maxLevel[i] - inventory[i][t];
            constraint delivery[t][i] <= maxLevel[i] * isDelivered[t][i];
        }
    }

    // Total inventory cost at the supplier
    costInventorySupplier <-
            holdingCostSupplier * sum[t in 0...horizonLength+1](inventorySupplier[t]);

    // Total inventory cost at customers
    costInventory <- sum[i in 0...nbCustomers][t in 0...horizonLength+1](
            holdingCost[i] * inventory[i][t]);

    // Total transportation cost
    costRoute <- sum[t in 0...horizonLength](costRoute[t]);

    // Objective: minimize the sum of all costs
    objective <- costInventorySupplier + costInventory + costRoute;
    minimize(objective);
}

/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 20;
}

/* Write the solution in a file with the following format :
 * - total distance run by the vehicle
 * - the nodes visited at each time step (omitting the start/end at the supplier) */
function output() {
    if (solFileName == nil) return;
    local outfile = io.openWrite(solFileName);

    outfile.println(costRoute.value);
    for [t in 0...horizonLength] {
        for [customer in route[t].value] {
            outfile.print(customer + 1, " ");
        }
        outfile.println();
    }
}

function readInputIrp() {
    inFile = io.openRead(inFileName);

    // General data
    firstLine = inFile.readln().trim().split();
    nbCustomers = firstLine[0].toInt()-1;
    horizonLength = firstLine[1].toInt();
    capacity = firstLine[2].toInt();

    // Supplier data
    secondLine = inFile.readln().trim().split();
    xCoordSupplier = secondLine[1].toDouble();
    yCoordSupplier = secondLine[2].toDouble();
    startLevelSupplier = secondLine[3].toInt();
    productionRate_supplier = secondLine[4].toInt();
    holdingCostSupplier = secondLine[5].toDouble();

    // Customers data
    for [i in 0...nbCustomers] {
        elements = inFile.readln().trim().split();
        xCoord[i] = elements[1].toDouble();
        yCoord[i] = elements[2].toDouble();
        startLevel[i] = elements[3].toInt();
        maxLevel[i] = elements[4].toInt();
        minLevel[i] = elements[5].toInt();
        demandRate[i] = elements[6].toInt();
        holdingCost[i] = elements[7].toDouble();
    }
}

function computeDistanceMatrix() {
    distance[i in 0...nbCustomers][j in 0...nbCustomers] =
            round(sqrt((xCoord[i] - xCoord[j]) * (xCoord[i] - xCoord[j])
            + (yCoord[i] - yCoord[j]) * (yCoord[i] - yCoord[j])));
    distanceSupplier[i in 0...nbCustomers] =
            round(sqrt((xCoord[i] - xCoordSupplier) * (xCoord[i] - xCoordSupplier)
            + (yCoord[i] - yCoordSupplier) * (yCoord[i] - yCoordSupplier)));
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python irp.py instances\abs3n10.dat
 
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python irp.py instances/abs3n10.dat
import hexaly.optimizer
import sys
import math


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


def main(instance_file, str_time_limit, sol_file):
    #
    # Read instance data
    #
    nb_customers, horizon_length, capacity, start_level_supplier, production_rate_supplier, \
        holding_cost_supplier, start_level, max_level, demand_rate, holding_cost, \
        dist_matrix_data, dist_supplier_data = read_input_irp(instance_file)

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

        # Quantity of product delivered at each discrete time instant of
        # the planning time horizon to each customer
        delivery = [[model.float(0, capacity) for i in range(nb_customers)]
                    for _ in range(horizon_length)]

        # Sequence of customers visited at each discrete time instant of
        # the planning time horizon
        route = [model.list(nb_customers) for t in range(horizon_length)]

        # Customers receive products only if they are visited
        is_delivered = [[model.contains(route[t], i) for i in range(nb_customers)]
                        for t in range(horizon_length)]

        # Create Hexaly arrays to be able to access them with an "at" operator
        dist_matrix = model.array(dist_matrix_data)
        dist_supplier = model.array(dist_supplier_data)

        dist_routes = [None for _ in range(horizon_length)]

        for t in range(horizon_length):
            sequence = route[t]
            c = model.count(sequence)

            # Distance traveled at instant t
            dist_lambda = model.lambda_function(
                lambda i:
                    model.at(
                        dist_matrix,
                        sequence[i - 1],
                        sequence[i]))
            dist_routes[t] = model.iif(
                c > 0,
                dist_supplier[sequence[0]]
                + model.sum(model.range(1, c), dist_lambda)
                + dist_supplier[sequence[c - 1]],
                0)

        # Stockout constraints at the supplier
        inventory_supplier = [None for _ in range(horizon_length + 1)]
        inventory_supplier[0] = start_level_supplier
        for t in range(1, horizon_length + 1):
            inventory_supplier[t] = inventory_supplier[t - 1] - model.sum(
                delivery[t - 1][i] for i in range(nb_customers)) + production_rate_supplier
            if t != horizon_length:
                model.constraint(
                    inventory_supplier[t] >= model.sum(delivery[t][i]
                                                       for i in range(nb_customers)))

        # Stockout constraints at the customers
        inventory = [[None for _ in range(horizon_length + 1)] for _ in range(nb_customers)]
        for i in range(nb_customers):
            inventory[i][0] = start_level[i]
            for t in range(1, horizon_length + 1):
                inventory[i][t] = inventory[i][t - 1] + delivery[t - 1][i] - demand_rate[i]
                model.constraint(inventory[i][t] >= 0)

        for t in range(horizon_length):
            # Capacity constraints
            model.constraint(
                model.sum((delivery[t][i]) for i in range(nb_customers)) <= capacity)

            # Maximum level constraints
            for i in range(nb_customers):
                model.constraint(delivery[t][i] <= max_level[i] - inventory[i][t])
                model.constraint(delivery[t][i] <= max_level[i] * is_delivered[t][i])

        # Total inventory cost at the supplier
        total_cost_inventory_supplier = holding_cost_supplier * model.sum(
            inventory_supplier[t] for t in range(horizon_length + 1))

        # Total inventory cost at customers
        total_cost_inventory = model.sum(model.sum(
            holding_cost[i] * inventory[i][t] for t in range(horizon_length + 1))
            for i in range(nb_customers))

        # Total transportation cost
        total_cost_route = model.sum(dist_routes[t] for t in range(horizon_length))

        # Objective: minimize the sum of all costs
        objective = total_cost_inventory_supplier + total_cost_inventory + total_cost_route
        model.minimize(objective)

        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 :
        # - total distance run by the vehicle
        # - the nodes visited at each time step (omitting the start/end at the supplier)
        #
        if len(sys.argv) >= 3:
            with open(sol_file, 'w') as f:
                f.write("%d\n" % (total_cost_route.value))
                for t in range(horizon_length):
                    for customer in route[t].value:
                        f.write("%d " % (customer + 1))
                    f.write("\n")


# The input files follow the "Archetti" format
def read_input_irp(filename):
    file_it = iter(read_elem(filename))

    nb_customers = int(next(file_it)) - 1
    horizon_length = int(next(file_it))
    capacity = int(next(file_it))

    x_coord = [None] * nb_customers
    y_coord = [None] * nb_customers
    start_level = [None] * nb_customers
    max_level = [None] * nb_customers
    min_level = [None] * nb_customers
    demand_rate = [None] * nb_customers
    holding_cost = [None] * nb_customers

    next(file_it)
    x_coord_supplier = float(next(file_it))
    y_coord_supplier = float(next(file_it))
    start_level_supplier = int(next(file_it))
    production_rate_supplier = int(next(file_it))
    holding_cost_supplier = float(next(file_it))
    for i in range(nb_customers):
        next(file_it)
        x_coord[i] = float(next(file_it))
        y_coord[i] = float(next(file_it))
        start_level[i] = int(next(file_it))
        max_level[i] = int(next(file_it))
        min_level[i] = int(next(file_it))
        demand_rate[i] = int(next(file_it))
        holding_cost[i] = float(next(file_it))

    distance_matrix = compute_distance_matrix(x_coord, y_coord)
    distance_supplier = compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord)

    return nb_customers, horizon_length, capacity, start_level_supplier, \
        production_rate_supplier, holding_cost_supplier, start_level, max_level, \
        demand_rate, holding_cost, distance_matrix, distance_supplier


# Compute the distance matrix
def compute_distance_matrix(x_coord, y_coord):
    nb_customers = len(x_coord)
    distance_matrix = [[None for i in range(nb_customers)] for j in range(nb_customers)]
    for i in range(nb_customers):
        distance_matrix[i][i] = 0
        for j in range(nb_customers):
            dist = compute_dist(x_coord[i], x_coord[j], y_coord[i], y_coord[j])
            distance_matrix[i][j] = dist
            distance_matrix[j][i] = dist
    return distance_matrix


# Compute the distances to the supplier
def compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord):
    nb_customers = len(x_coord)
    distance_supplier = [None] * nb_customers
    for i in range(nb_customers):
        dist = compute_dist(x_coord_supplier, x_coord[i], y_coord_supplier, y_coord[i])
        distance_supplier[i] = dist
    return distance_supplier


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


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

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

    main(instance_file, str_time_limit, sol_file)
Compilation / Execution (Windows)
cl /EHsc irp.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
irp instances\abs3n10.dat
 
Compilation / Execution (Linux)
g++ irp.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o irp
./irp instances/abs3n10.dat
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

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

    // Number of customers
    int nbCustomers;

    // Horizon length
    int horizonLength;

    // Capacity
    int capacity;

    // Start level at the supplier
    hxint startLevelSupplier;

    // Production rate of the supplier
    int productionRateSupplier;

    // Holding costs of the supplier
    double holdingCostSupplier;

    // Start level of the customers
    vector<hxint> startLevel;

    // Max level of the customers
    vector<int> maxLevel;

    // Demand rate of the customers
    vector<int> demandRate;

    // Holding costs of the customers
    vector<double> holdingCost;

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

    // Distance to depot
    vector<int> distSupplierData;

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

    // Decision variables
    vector<HxExpression> route;

    // Are the customers receiving products
    vector<vector<HxExpression>> isDelivered;

    // Total inventory cost at the supplier
    HxExpression totalCostInventorySupplier;

    // Total inventory cost at customers
    HxExpression totalCostInventory;

    // Total transportation cost
    HxExpression totalCostRoute;

    // Objective
    HxExpression objective;

    Irp() {}

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

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

        // Quantity of product delivered at each discrete time instant of
        // the planning time horizon to each customer
        delivery.resize(horizonLength);
        for (int t = 0; t < horizonLength; ++t) {
            delivery[t].resize(nbCustomers);
            for (int i = 0; i < nbCustomers; ++i) {
                delivery[t][i] = model.floatVar(0, capacity);
            }
        }

        // Sequence of customers visited at each discrete time instant of
        // the planning time horizon
        route.resize(horizonLength);
        for (int t = 0; t < horizonLength; ++t) {
            route[t] = model.listVar(nbCustomers);
        }

        // Create distances as arrays to be able to access them with an "at" operator
        HxExpression distMatrix = model.array();
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrix.addOperand(model.array(distMatrixData[i].begin(), distMatrixData[i].end()));
        }
        HxExpression distSupplier = model.array(distSupplierData.begin(), distSupplierData.end());

        isDelivered.resize(horizonLength);
        vector<HxExpression> distRoutes(horizonLength);

        for (int t = 0; t < horizonLength; ++t) {
            HxExpression sequence = route[t];
            HxExpression c = model.count(sequence);
            isDelivered[t].resize(nbCustomers);

            // Customers receive products only if they are visited
            for (int i = 0; i < nbCustomers; ++i) {
                isDelivered[t][i] = model.contains(sequence, i);
            }

            // Distance traveled at instant t
            HxExpression distLambda = model.createLambdaFunction(
                [&](HxExpression i) { return model.at(distMatrix, sequence[i - 1], sequence[i]); });
            distRoutes[t] = model.iif(c > 0,
                                     distSupplier[sequence[0]] + model.sum(model.range(1, c), distLambda) +
                                         distSupplier[sequence[c - 1]],
                                     0);
        }

        // Stockout constraints at the supplier
        vector<HxExpression> inventorySupplier(horizonLength + 1);
        inventorySupplier[0] = model.createConstant(startLevelSupplier);
        for (int t = 1; t < horizonLength + 1; ++t) {
            inventorySupplier[t] = inventorySupplier[t - 1] -
                                   model.sum(delivery[t - 1].begin(), delivery[t - 1].end()) + productionRateSupplier;
            if (t != horizonLength) {
                model.constraint(inventorySupplier[t] >= model.sum(delivery[t].begin(), delivery[t].end()));
            }
        }

        // Stockout constraints at the customers
        vector<vector<HxExpression>> inventory(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            inventory[i].resize(horizonLength + 1);
            inventory[i][0] = model.createConstant(startLevel[i]);
            for (int t = 1; t < horizonLength + 1; ++t) {
                inventory[i][t] = inventory[i][t - 1] + delivery[t - 1][i] - demandRate[i];
                model.constraint(inventory[i][t] >= 0);
            }
        }

        for (int t = 0; t < horizonLength; ++t) {
            // Capacity constraints
            model.constraint(model.sum(delivery[t].begin(), delivery[t].end()) <= capacity);

            // Maximum level constraints
            for (int i = 0; i < nbCustomers; ++i) {
                model.constraint(delivery[t][i] <= maxLevel[i] - inventory[i][t]);
                model.constraint(delivery[t][i] <= maxLevel[i] * isDelivered[t][i]);
            }
        }

        // Total inventory cost at the supplier
        totalCostInventorySupplier =
            holdingCostSupplier * model.sum(inventorySupplier.begin(), inventorySupplier.end());

        // Total inventory cost at customers
        vector<HxExpression> costInventoryCustomer(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            costInventoryCustomer[i] = holdingCost[i] * model.sum(inventory[i].begin(), inventory[i].end());
        }
        totalCostInventory = model.sum(costInventoryCustomer.begin(), costInventoryCustomer.end());

        // Total transportation cost
        totalCostRoute = model.sum(distRoutes.begin(), distRoutes.end());

        // Objective: minimize the sum of all costs
        objective = totalCostInventorySupplier + totalCostInventory + totalCostRoute;
        model.minimize(objective);

        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     * - total distance run by the vehicle
     * - the nodes visited at each time step (omitting the start/end at the supplier) */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << totalCostRoute.getValue() << endl;
        for (int t = 0; t < horizonLength; ++t) {
            HxCollection routeCollection = route[t].getCollectionValue();
            for (int i = 0; i < routeCollection.count(); ++i) {
                outfile << routeCollection[i] + 1 << " ";
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the "Archetti" format
    void readInputIrp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }

        infile >> nbCustomers;
        nbCustomers -= 1;
        infile >> horizonLength;
        infile >> capacity;
        int idSupplier;
        double xCoordSupplier, yCoordSupplier;
        vector<int> id;
        vector<double> xCoord, yCoord;
        infile >> idSupplier;
        infile >> xCoordSupplier;
        infile >> yCoordSupplier;
        infile >> startLevelSupplier;
        infile >> productionRateSupplier;
        infile >> holdingCostSupplier;
        vector<int> minLevel;
        id.resize(nbCustomers);
        xCoord.resize(nbCustomers);
        yCoord.resize(nbCustomers);
        startLevel.resize(nbCustomers);
        maxLevel.resize(nbCustomers);
        minLevel.resize(nbCustomers);
        demandRate.resize(nbCustomers);
        holdingCost.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            infile >> id[i];
            infile >> xCoord[i];
            infile >> yCoord[i];
            infile >> startLevel[i];
            infile >> maxLevel[i];
            infile >> minLevel[i];
            infile >> demandRate[i];
            infile >> holdingCost[i];
        }

        printf("%lf", holdingCost[3]);

        computeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);

        infile.close();
    }

    // Compute the distance matrix
    void computeDistanceMatrix(double xCoordSupplier, double yCoordSupplier, const vector<double>& xCoord,
                               const vector<double>& yCoord) {
        distMatrixData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i].resize(nbCustomers);
        }
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                int distance = computeDist(xCoord[i], xCoord[j], yCoord[i], yCoord[j]);
                distMatrixData[i][j] = distance;
                distMatrixData[j][i] = distance;
            }
        }

        distSupplierData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distSupplierData[i] = computeDist(xCoordSupplier, xCoord[i], yCoordSupplier, yCoord[i]);
        }
    }

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

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: irp 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 {
        Irp model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc Irp.cs /reference:Hexaly.NET.dll
Irp instances\abs3n10.dat
using System;
using System.IO;
using System.Globalization;
using Hexaly.Optimizer;

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

    // Number of customers
    int nbCustomers;

    // Horizon length
    int horizonLength;

    // Capacity
    int capacity;

    // Start level at the supplier
    int startLevelSupplier;

    // Production rate of the supplier
    int productionRateSupplier;

    // Holding costs of the supplier
    double holdingCostSupplier;

    // Start level of the customers
    int[] startLevel;

    // Max level of the customers
    int[] maxLevel;

    // Demand rate of the customers
    int[] demandRate;

    // Holding costs of the customers
    double[] holdingCost;

    // Distance matrix between customers
    long[][] distMatrixData;

    // Distances between customers and supplier
    long[] distSupplierData;

    // Decision variables
    private HxExpression[][] delivery;

    // Decision variables
    private HxExpression[] route;

    // Are the customers receiving products
    private HxExpression[][] isDelivered;

    // Distance traveled by each truck
    HxExpression[] distRoutes;

    // Inventory at the supplier
    HxExpression[] inventorySupplier;

    // Inventory at the customers
    HxExpression[][] inventory;

    // Total inventory cost at the supplier
    HxExpression totalCostInventorySupplier;

    // Inventory cost at a customer
    HxExpression[] costInventoryCustomer;

    // Total inventory cost at customers
    HxExpression totalCostInventory;

    // Total transportation cost
    HxExpression totalCostRoute;

    // Objective
    HxExpression objective;

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

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

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

        delivery = new HxExpression[horizonLength][];
        route = new HxExpression[horizonLength];
        isDelivered = new HxExpression[horizonLength][];
        distRoutes = new HxExpression[horizonLength];

        // Quantity of product delivered at each discrete time instant of
        // the planning time horizon to each customer
        for (int t = 0; t < horizonLength; ++t)
        {
            delivery[t] = new HxExpression[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
                delivery[t][i] = model.Float(0, capacity);
        }

        // Sequence of customers visited at each discrete time instant of
        // the planning time horizon
        for (int t = 0; t < horizonLength; ++t)
            route[t] = model.List(nbCustomers);

        // Create distances as arrays to be able to access them with an "at" operator
        HxExpression distSupplier = model.Array(distSupplierData);
        HxExpression distMatrix = model.Array(distMatrixData);

        for (int t = 0; t < horizonLength; ++t)
        {
            HxExpression sequence = route[t];
            HxExpression c = model.Count(sequence);

            isDelivered[t] = new HxExpression[nbCustomers];
            // Customers receive products only if they are visited
            for (int i = 0; i < nbCustomers; ++i)
                isDelivered[t][i] = model.Contains(sequence, i);

            // Distance traveled at instant t
            HxExpression distLambda = model.LambdaFunction(
                i => distMatrix[sequence[i - 1], sequence[i]]
            );
            distRoutes[t] = model.If(
                c > 0,
                distSupplier[sequence[0]]
                    + model.Sum(model.Range(1, c), distLambda)
                    + distSupplier[sequence[c - 1]],
                0
            );
        }

        // Stockout constraints at the supplier
        inventorySupplier = new HxExpression[horizonLength + 1];
        inventorySupplier[0] = model.CreateConstant(startLevelSupplier);
        for (int t = 1; t < horizonLength + 1; ++t)
        {
            inventorySupplier[t] =
                inventorySupplier[t - 1] - model.Sum(delivery[t - 1]) + productionRateSupplier;
            if (t != horizonLength)
                model.Constraint(inventorySupplier[t] >= model.Sum(delivery[t]));
        }

        // Stockout constraints at the customers
        inventory = new HxExpression[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
        {
            inventory[i] = new HxExpression[horizonLength + 1];
            inventory[i][0] = model.CreateConstant(startLevel[i]);
            for (int t = 1; t < horizonLength + 1; ++t)
            {
                inventory[i][t] = inventory[i][t - 1] + delivery[t - 1][i] - demandRate[i];
                model.Constraint(inventory[i][t] >= 0);
            }
        }

        for (int t = 0; t < horizonLength; ++t)
        {
            // Capacity constraints
            model.Constraint(model.Sum(delivery[t]) <= capacity);

            // Maximum level constraints
            for (int i = 0; i < nbCustomers; ++i)
            {
                model.Constraint(delivery[t][i] <= maxLevel[i] - inventory[i][t]);
                model.Constraint(delivery[t][i] <= maxLevel[i] * isDelivered[t][i]);
            }
        }

        // Total inventory cost at the supplier
        totalCostInventorySupplier = holdingCostSupplier * model.Sum(inventorySupplier);

        // Total inventory cost at customers
        costInventoryCustomer = new HxExpression[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            costInventoryCustomer[i] = holdingCost[i] * model.Sum(inventory[i]);
        totalCostInventory = model.Sum(costInventoryCustomer);

        // Total transportation cost
        totalCostRoute = model.Sum(distRoutes);

        // Objective: minimize the sum of all costs
        objective = model.Sum(
            model.Sum(totalCostInventorySupplier, totalCostInventory),
            totalCostRoute
        );
        model.Minimize(objective);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     * - total distance run by the vehicle
     * - the nodes visited at each time step (omitting the start/end at the supplier) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalCostRoute.GetValue());
            for (int t = 0; t < horizonLength; ++t)
            {
                HxCollection routeCollection = route[t].GetCollectionValue();
                for (int i = 0; i < routeCollection.Count(); ++i)
                    output.Write((routeCollection[i] + 1) + " ");
                output.WriteLine();
            }
        }
    }

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

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

    // The input files follow the "Archetti" format
    private void ReadInputIrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            splitted = ReadNextLine(input);
            nbCustomers = int.Parse(splitted[0]) - 1;
            horizonLength = int.Parse(splitted[1]);
            capacity = int.Parse(splitted[2]);
            splitted = ReadNextLine(input);
            Console.WriteLine(1);
            Console.WriteLine(splitted[1]);
            double xCoordSupplier = double.Parse(splitted[1], CultureInfo.InvariantCulture);
            Console.WriteLine(2);
            double yCoordSupplier = double.Parse(splitted[2], CultureInfo.InvariantCulture);
            Console.WriteLine(3);
            startLevelSupplier = int.Parse(splitted[3]);
            Console.WriteLine(4);
            productionRateSupplier = int.Parse(splitted[4]);
            holdingCostSupplier = double.Parse(splitted[5], CultureInfo.InvariantCulture);
            double[] xCoord = new double[nbCustomers];
            double[] yCoord = new double[nbCustomers];

            startLevel = new int[nbCustomers];
            maxLevel = new int[nbCustomers];
            demandRate = new int[nbCustomers];
            holdingCost = new double[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i)
            {
                splitted = ReadNextLine(input);
                xCoord[i] = double.Parse(splitted[1], CultureInfo.InvariantCulture);
                yCoord[i] = double.Parse(splitted[2], CultureInfo.InvariantCulture);
                startLevel[i] = int.Parse(splitted[3]);
                maxLevel[i] = int.Parse(splitted[4]);
                demandRate[i] = int.Parse(splitted[6]);
                holdingCost[i] = double.Parse(splitted[7], CultureInfo.InvariantCulture);
            }

            ComputeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);
        }
    }

    // Compute the distance matrix
    private void ComputeDistanceMatrix(
        double xCoordSupplier,
        double yCoordSupplier,
        double[] xCoord,
        double[] yCoord
    )
    {
        distMatrixData = new long[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new long[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                long dist = ComputeDist(xCoord[i], xCoord[j], yCoord[i], yCoord[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distSupplierData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            distSupplierData[i] = ComputeDist(xCoordSupplier, xCoord[i], yCoordSupplier, yCoord[i]);
    }

    private long ComputeDist(double xi, double xj, double yi, double yj)
    {
        return Convert.ToInt64(Math.Round(Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2))));
    }
}
Compilation / Execution (Windows)
javac Irp.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. Irp instances\abs3n10.dat
Compilation / Execution (Linux)
javac Irp.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. Irp instances/abs3n10.dat
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

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

    // Number of customers
    private int nbCustomers;

    // Horizon length
    private int horizonLength;

    // Capacity
    private int capacity;

    // Start level at the supplier
    private int startLevelSupplier;

    // Production rate of the supplier
    private int productionRateSupplier;

    // Holding costs of the supplier
    private double holdingCostSupplier;

    // Start level of the customers
    private int[] startLevel;

    // Max level of the customers
    private int[] maxLevel;

    // Demand rate of the customers
    private int[] demandRate;

    // Holding costs of the customers
    private double[] holdingCost;

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

    // Distances between customers and supplier
    private long[] distSupplierData;

    // Decision variables
    private HxExpression[][] delivery;

    // Decision variables
    private HxExpression[] route;

    // Are the customers receiving products
    private HxExpression[][] isDelivered;

    // Distance traveled by each truck
    private HxExpression[] distRoutes;

    // Inventory at the supplier
    private HxExpression[] inventorySupplier;

    // Inventory at the customers
    private HxExpression[][] inventory;

    // Total inventory cost at the supplier
    private HxExpression totalCostInventorySupplier;

    // Inventory cost at a customer
    private HxExpression[] costInventoryCustomer;

    // Total inventory cost at customers
    private HxExpression totalCostInventory;

    // Total transportation cost
    private HxExpression totalCostRoute;

    // Objective
    private HxExpression objective;

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

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

        delivery = new HxExpression[horizonLength][nbCustomers];
        route = new HxExpression[horizonLength];
        isDelivered = new HxExpression[horizonLength][nbCustomers];
        distRoutes = new HxExpression[horizonLength];

        // Quantity of product delivered at each discrete time instant of
        // the planning time horizon to each customer
        for (int t = 0; t < horizonLength; ++t) {
            for (int i = 0; i < nbCustomers; ++i) {
                delivery[t][i] = model.floatVar(0, capacity);
            }
        }

        // Sequence of customers visited at each discrete time instant of
        // the planning time horizon
        for (int t = 0; t < horizonLength; ++t) {
            route[t] = model.listVar(nbCustomers);
        }

        // Create distances as arrays to be able to access them with an "at" operator
        HxExpression distSupplier = model.array(distSupplierData);
        HxExpression distMatrix = model.array(distMatrixData);

        for (int t = 0; t < horizonLength; ++t) {
            HxExpression sequence = route[t];
            HxExpression c = model.count(sequence);

            // Customers receive products only if they are visited
            for (int i = 0; i < nbCustomers; ++i) {
                isDelivered[t][i] = model.contains(sequence, i);
            }

            // Distance traveled at instant t
            HxExpression distLambda = model
                .lambdaFunction(i -> model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)));
            distRoutes[t] = model.iif(model.gt(c, 0),
                model.sum(
                    model.sum(model.at(distSupplier, model.at(sequence, 0)), model.sum(model.range(1, c), distLambda)),
                    model.at(distSupplier, model.at(sequence, model.sub(c, 1)))),
                0);
        }

        // Stockout constraints at the supplier
        inventorySupplier = new HxExpression[horizonLength + 1];
        inventorySupplier[0] = model.createConstant(startLevelSupplier);
        for (int t = 1; t < horizonLength + 1; ++t) {
            inventorySupplier[t] = model.sum(model.sub(inventorySupplier[t - 1], model.sum(delivery[t - 1])),
                productionRateSupplier);
            if (t != horizonLength) {
                model.constraint(model.geq(inventorySupplier[t], model.sum(delivery[t])));
            }
        }

        // Stockout constraints at the customers
        inventory = new HxExpression[nbCustomers][horizonLength + 1];
        for (int i = 0; i < nbCustomers; ++i) {
            inventory[i][0] = model.createConstant(startLevel[i]);
            for (int t = 1; t < horizonLength + 1; ++t) {
                inventory[i][t] = model.sub(model.sum(inventory[i][t - 1], delivery[t - 1][i]), demandRate[i]);
                model.constraint(model.geq(inventory[i][t], 0));
            }
        }

        for (int t = 0; t < horizonLength; ++t) {
            // Capacity constraints
            model.constraint(model.leq(model.sum(delivery[t]), capacity));

            // Maximum level constraints
            for (int i = 0; i < nbCustomers; ++i) {
                model.constraint(model.leq(delivery[t][i], model.sub(maxLevel[i], inventory[i][t])));
                model.constraint(model.leq(delivery[t][i], model.prod(maxLevel[i], isDelivered[t][i])));
            }
        }

        // Total inventory cost at the supplier
        totalCostInventorySupplier = model.prod(holdingCostSupplier, model.sum(inventorySupplier));

        // Total inventory cost at customers
        costInventoryCustomer = new HxExpression[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            costInventoryCustomer[i] = model.prod(holdingCost[i], model.sum(inventory[i]));
        }
        totalCostInventory = model.sum(costInventoryCustomer);

        // Total transportation cost
        totalCostRoute = model.sum(distRoutes);

        // Objective: minimize the sum of all costs
        objective = model.sum(model.sum(totalCostInventorySupplier, totalCostInventory), totalCostRoute);
        model.minimize(objective);

        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     * - total distance run by the vehicle
     * - the nodes visited at each time step (omitting the start/end at the supplier) */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println(totalCostRoute.getValue());
            for (int t = 0; t < horizonLength; ++t) {
                HxCollection routeCollection = route[t].getCollectionValue();
                for (int i = 0; i < routeCollection.count(); ++i) {
                    output.print((routeCollection.get(i) + 1) + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "Archetti" format
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            String[] splitted;
            splitted = input.nextLine().split("\\s+");
            nbCustomers = Integer.parseInt(splitted[0]) - 1;
            horizonLength = Integer.parseInt(splitted[1]);
            capacity = Integer.parseInt(splitted[2]);
            splitted = input.nextLine().split("\\s+");
            double xCoordSupplier = Double.parseDouble(splitted[2]);
            double yCoordSupplier = Double.parseDouble(splitted[3]);
            startLevelSupplier = Integer.parseInt(splitted[4]);
            productionRateSupplier = Integer.parseInt(splitted[5]);
            holdingCostSupplier = Double.parseDouble(splitted[6]);
            double[] xCoord = new double[nbCustomers];
            double[] yCoord = new double[nbCustomers];
            startLevel = new int[nbCustomers];
            maxLevel = new int[nbCustomers];
            demandRate = new int[nbCustomers];
            holdingCost = new double[nbCustomers];
            for (int i = 0; i < nbCustomers; ++i) {
                splitted = input.nextLine().split("\\s+");
                xCoord[i] = Double.parseDouble(splitted[2]);
                yCoord[i] = Double.parseDouble(splitted[3]);
                startLevel[i] = Integer.parseInt(splitted[4]);
                maxLevel[i] = Integer.parseInt(splitted[5]);
                demandRate[i] = Integer.parseInt(splitted[7]);
                holdingCost[i] = Double.parseDouble(splitted[8]);
            }

            computeDistanceMatrix(xCoordSupplier, yCoordSupplier, xCoord, yCoord);
        }
    }

    // Compute the distance matrix
    private void computeDistanceMatrix(double xCoordSupplier, double yCoordSupplier, double[] xCoord, double[] yCoord) {
        distMatrixData = new long[nbCustomers][nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                long dist = computeDist(xCoord[i], xCoord[j], yCoord[i], yCoord[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distSupplierData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distSupplierData[i] = computeDist(xCoordSupplier, xCoord[i], yCoordSupplier, yCoord[i]);
        }
    }

    private long computeDist(double xi, double xj, double yi, double yj) {
        return Math.round(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 Irp 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";

            Irp model = new Irp(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}

Results

Hexaly Optimizer improves the literature’s best known solutions on 36% of IRP instances from the 2016 ROADEF/EURO Challenge. Our Inventory Routing Problem (IRP) benchmark page illustrates Hexaly Optimizer’s performance on this challenging problem.

Similar examples

Discover the ease of use and performance of Hexaly