Pickup and Delivery Problem with Time Windows (PDPTW)
Problem
In the Pickup and Delivery Problem with Time Windows (PDPTW), a fleet of delivery vehicles with uniform capacity must collect and deliver items according to the customers’ demands and opening hours. The customers are paired: in each pair, one customer corresponds to a pickup point (with positive demand), and the other corresponds to a delivery point (with negative demand). Each customer must be served by exactly one vehicle. Furthermore, paired customers must be served by the same truck and each pickup must be treated before its associated delivery. The vehicles start and end their routes at a common depot and the load they carry must not exceed their capacity at any point in the tour. The objectives consist in minimizing the fleet size and the total traveled distance.
Principles learned
- Add list decision variables to model the trucks’ sequences of customers
- Use the ‘find’ operator to ensure paired customers are served by the same truck
- Define an array using a recursive lambda function to compute the trucks’ load over time and the customers’ visiting times
Data
The instances provided come from the Li & Lim benchmark. The format of the data files is as follows:
- The first line gives the number of vehicles, the capacity, and the speed (not used)
- From the second line, for each customer (starting with the depot):
- The index of the customer
- The x coordinate
- The y coordinate
- The demand
- The earliest arrival
- The latest arrival
- The service time
- The index of the corresponding pickup order (0 if the order is a delivery)
- The index of the corresponding delivery order (0 if the order is a pickup)
Program
The Hexaly model for the Pickup and Delivery Problem with Time Windows (PDPTW) is an extension of the CVRPTW model. We refer the reader to this model for the routing and time-window aspects of the problem.
Paired customers must be visited by the same truck. Using the ‘find’ operator, we retrieve the two list variables containing the pickup point and delivery point of each pair respectively. We can then constrain these two lists to be the same. Furthermore, the goods must be picked up before they are delivered. Using the ‘indexOf’ operator, we ensure that the pickup point is placed before the delivery point in the list.
The load each truck carries varies along the tour: it increases upon pickups and then decreases upon deliveries. We compute the trucks’ load over time using a recursive array: the load after visiting a customer is equal to the load after visiting the previous customer plus the current customer’s demand (positive in case of a pickup, negative in case of a delivery). We then use a variadic ‘and’ operator to ensure the capacity is respected at any point in the tour.
Finally, the objectives are the same as for the CVRPTW: we minimize the total lateness, the number of trucks used, and the total traveled distance.
- Execution
-
hexaly pdptw.hxm inFileName=instances/lc101.txt [hxTimeLimit=] [solFileName=]
use io;
/* Read instance data. The input files follow the "Li & Lim" format*/
function input() {
usage = "Usage: hexaly pdptw.hxm "
+ "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";
if (inFileName == nil) throw usage;
readInputPdptw();
// Compute distance matrix
computeDistanceMatrix();
if (nbMaxTrucks != nil) {
nbTrucks = nbMaxTrucks;
}
}
/* Declare the optimization model */
function model() {
customersSequences[k in 0...nbTrucks] <- list(nbCustomers);
// All customers must be visited by exactly one truck
constraint partition[k in 0...nbTrucks](customersSequences[k]);
//Pickups and deliveries
for[i in 0...nbCustomers : pickupIndex[i] == -1] {
pickupListIndex <- find(customersSequences, i);
deliveryListIndex <- find(customersSequences, deliveryIndex[i]);
constraint pickupListIndex == deliveryListIndex;
pickupList <- customersSequences[pickupListIndex];
deliveryList <- customersSequences[deliveryListIndex];
constraint indexOf(pickupList, i) < indexOf(deliveryList, deliveryIndex[i]);
}
for [k in 0...nbTrucks] {
local sequence <- customersSequences[k];
local c <- count(sequence);
// A truck is used if it visits at least one customer
truckUsed[k] <- c > 0;
// The quantity needed in each route must not exceed the truck capacity
// at any point in the sequence
routeQuantity[k] <- array(0...c, (i, prev) => prev + demands[sequence[i]], 0);
constraint and(0...c, i => routeQuantity[k][i] <= truckCapacity);
endTime[k] <- array(0...c, (i, prev) => max(earliestStart[sequence[i]],
i == 0 ? distanceDepot[sequence[0]] :
prev + distanceMatrix[sequence[i - 1]][sequence[i]])
+ serviceTime[sequence[i]], 0);
homeLateness[k] <- truckUsed[k] ?
max(0, endTime[k][c - 1] + distanceDepot[sequence[c - 1]] - maxHorizon) :
0;
// Distance traveled by truck k
routeDistances[k] <- sum(1...c,
i => distanceMatrix[sequence[i - 1]][sequence[i]]) + (truckUsed[k] ?
(distanceDepot[sequence[0]] + distanceDepot[sequence[c - 1]]) :
0);
lateness[k] <- homeLateness[k] + sum(0...c,
i => max(0, endTime[k][i] - latestEnd[sequence[i]]));
}
// Total lateness, must be 0 for a solution to be valid
totalLateness <- sum[k in 0...nbTrucks](lateness[k]);
nbTrucksUsed <- sum[k in 0...nbTrucks](truckUsed[k]);
// Total distance traveled
totalDistance <- round(100 * sum[k in 0...nbTrucks](routeDistances[k])) / 100;
minimize totalLateness;
if (nbMaxTrucks == nil) minimize nbTrucksUsed;
minimize totalDistance;
}
/* Parametrize the solver */
function param() {
if (hxTimeLimit == nil) hxTimeLimit = 20;
}
/* Write the solution in a file with the following format :
* - number of trucks used and total distance
* - for each truck the customers visited (omitting the start/end at the depot) */
function output() {
if (solFileName == nil) return;
local outfile = io.openWrite(solFileName);
outfile.println(nbTrucksUsed.value, " ", totalDistance.value);
for [k in 0...nbTrucks] {
if (truckUsed[k].value != 1) continue;
for [customer in customersSequences[k].value]
outfile.print(customer + 1, " ");
outfile.println();
}
}
function readInputPdptw() {
local inFile = io.openRead(inFileName);
// Truck related data
nbTrucks = inFile.readInt();
truckCapacity = inFile.readInt();
speed = inFile.readInt(); // not used
// Depot data
local line = inFile.readln().split(" ");
depotIndex = line[0].toInt();
depotX = line[1].toInt();
depotY = line[2].toInt();
maxHorizon = line[5].toInt();
// Customers data
i = 0;
while (!inFile.eof()) {
inLine = inFile.readln();
line = inLine.split(" ");
if (count(line) == 0) break;
if (count(line) != 9) throw "Wrong file format";
customerIndex[i] = line[0].toInt();
customerX[i] = line[1].toInt();
customerY[i] = line[2].toInt();
demands[i] = line[3].toInt();
serviceTime[i] = line[6].toInt();
earliestStart[i] = line[4].toInt();
// in input files due date is meant as latest start time
latestEnd[i] = line[5].toInt() + serviceTime[i];
pickupIndex[i] = line[7].toInt() - 1;
deliveryIndex[i] = line[8].toInt() - 1;
i = i + 1;
}
nbCustomers = i;
inFile.close();
}
function computeDistanceMatrix() {
for [i in 0...nbCustomers] {
distanceMatrix[i][i] = 0;
for [j in i+1...nbCustomers] {
local localDistance = computeDist(i, j);
distanceMatrix[j][i] = localDistance;
distanceMatrix[i][j] = localDistance;
}
}
for [i in 0...nbCustomers] {
local localDistance = computeDepotDist(i);
distanceDepot[i] = localDistance;
}
}
function computeDist(i, j) {
local x1 = customerX[i];
local x2 = customerX[j];
local y1 = customerY[i];
local y2 = customerY[j];
return computeDistance(x1, x2, y1, y2);
}
function computeDepotDist(i) {
local x1 = customerX[i];
local xd = depotX;
local y1 = customerY[i];
local yd = depotY;
return computeDistance(x1, xd, y1, yd);
}
function computeDistance(x1, x2, y1, y2) {
return sqrt(pow((x1 - x2), 2) + pow((y1 - y2), 2));
}
- Execution (Windows)
-
set PYTHONPATH=%HX_HOME%\bin\pythonpython pdptw.py instances\lc101.txt
- Execution (Linux)
-
export PYTHONPATH=/opt/hexaly_13_0/bin/pythonpython pdptw.py instances/lc101.txt
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_trucks, truck_capacity, dist_matrix_data, dist_depot_data, \
demands_data, service_time_data, earliest_start_data, latest_end_data, \
pick_up_index, delivery_index, max_horizon = read_input_pdptw(instance_file)
with hexaly.optimizer.HexalyOptimizer() as optimizer:
#
# Declare the optimization model
#
model = optimizer.model
# Sequence of customers visited by each truck
customers_sequences = [model.list(nb_customers) for k in range(nb_trucks)]
# All customers must be visited by exactly one truck
model.constraint(model.partition(customers_sequences))
# /Create Hexaly arrays to be able to access them with "at" operators
demands = model.array(demands_data)
earliest = model.array(earliest_start_data)
latest = model.array(latest_end_data)
service_time = model.array(service_time_data)
dist_matrix = model.array(dist_matrix_data)
dist_depot = model.array(dist_depot_data)
dist_routes = [None] * nb_trucks
end_time = [None] * nb_trucks
home_lateness = [None] * nb_trucks
lateness = [None] * nb_trucks
# A truck is used if it visits at least one customer
trucks_used = [(model.count(customers_sequences[k]) > 0) for k in range(nb_trucks)]
nb_trucks_used = model.sum(trucks_used)
# Pickups and deliveries
customers_sequences_array = model.array(customers_sequences)
for i in range(nb_customers):
if pick_up_index[i] == -1:
pick_up_list_index = model.find(customers_sequences_array, i)
delivery_list_index = model.find(customers_sequences_array, delivery_index[i])
model.constraint(pick_up_list_index == delivery_list_index)
pick_up_list = model.at(customers_sequences_array, pick_up_list_index)
delivery_list = model.at(customers_sequences_array, delivery_list_index)
model.constraint(model.index(pick_up_list, i) < model.index(delivery_list, delivery_index[i]))
for k in range(nb_trucks):
sequence = customers_sequences[k]
c = model.count(sequence)
# The quantity needed in each route must not exceed the truck capacity at any
# point in the sequence
demand_lambda = model.lambda_function(
lambda i, prev: prev + demands[sequence[i]])
route_quantity = model.array(model.range(0, c), demand_lambda, 0)
quantity_lambda = model.lambda_function(
lambda i: route_quantity[i] <= truck_capacity)
model.constraint(model.and_(model.range(0, c), quantity_lambda))
# Distance traveled by each truck
dist_lambda = model.lambda_function(
lambda i: model.at(dist_matrix, sequence[i - 1], sequence[i]))
dist_routes[k] = model.sum(model.range(1, c), dist_lambda) \
+ model.iif(c > 0, dist_depot[sequence[0]] + dist_depot[sequence[c - 1]], 0)
# End of each visit
end_lambda = model.lambda_function(
lambda i, prev:
model.max(
earliest[sequence[i]],
model.iif(
i == 0,
dist_depot[sequence[0]],
prev + model.at(dist_matrix, sequence[i - 1], sequence[i])))
+ service_time[sequence[i]])
end_time[k] = model.array(model.range(0, c), end_lambda, 0)
# Arriving home after max_horizon
home_lateness[k] = model.iif(
trucks_used[k],
model.max(
0,
end_time[k][c - 1] + dist_depot[sequence[c - 1]] - max_horizon),
0)
# Completing visit after latest_end
late_selector = model.lambda_function(
lambda i: model.max(0, end_time[k][i] - latest[sequence[i]]))
lateness[k] = home_lateness[k] + model.sum(model.range(0, c), late_selector)
# Total lateness (must be 0 for the solution to be valid)
total_lateness = model.sum(lateness)
# Total distance traveled
total_distance = model.div(model.round(100 * model.sum(dist_routes)), 100)
# Objective: minimize the number of trucks used, then minimize the distance traveled
model.minimize(total_lateness)
model.minimize(nb_trucks_used)
model.minimize(total_distance)
model.close()
# Parameterize the optimizer
optimizer.param.time_limit = int(str_time_limit)
optimizer.solve()
#
# Write the solution in a file with the following format:
# - number of trucks used and total distance
# - for each truck the customers visited (omitting the start/end at the depot)
#
if sol_file is not None:
with open(sol_file, 'w') as f:
f.write("%d %.2f\n" % (nb_trucks_used.value, total_distance.value))
for k in range(nb_trucks):
if trucks_used[k].value != 1:
continue
# Values in sequence are in 0...nbCustomers. +2 is to put it back in
# 2...nbCustomers+2 as in the data files (1 being the depot)
for customer in customers_sequences[k].value:
f.write("%d " % (customer + 1))
f.write("\n")
# The input files follow the "Li & Lim" format
def read_input_pdptw(filename):
file_it = iter(read_elem(filename))
nb_trucks = int(next(file_it))
truck_capacity = int(next(file_it))
next(file_it)
next(file_it)
depot_x = int(next(file_it))
depot_y = int(next(file_it))
for i in range(2):
next(file_it)
max_horizon = int(next(file_it))
for i in range(3):
next(file_it)
customers_x = []
customers_y = []
demands = []
earliest_start = []
latest_end = []
service_time = []
pick_up_index = []
delivery_index = []
while True:
val = next(file_it, None)
if val is None:
break
i = int(val) - 1
customers_x.append(int(next(file_it)))
customers_y.append(int(next(file_it)))
demands.append(int(next(file_it)))
ready = int(next(file_it))
due = int(next(file_it))
stime = int(next(file_it))
pick = int(next(file_it))
delivery = int(next(file_it))
earliest_start.append(ready)
# in input files due date is meant as latest start time
latest_end.append(due + stime)
service_time.append(stime)
pick_up_index.append(pick - 1)
delivery_index.append(delivery - 1)
nb_customers = i + 1
distance_matrix = compute_distance_matrix(customers_x, customers_y)
distance_depots = compute_distance_depots(depot_x, depot_y, customers_x, customers_y)
return nb_customers, nb_trucks, truck_capacity, distance_matrix, distance_depots, \
demands, service_time, earliest_start, latest_end, pick_up_index, \
delivery_index, max_horizon
# Compute the distance matrix
def compute_distance_matrix(customers_x, customers_y):
nb_customers = len(customers_x)
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(customers_x[i], customers_x[j],
customers_y[i], customers_y[j])
distance_matrix[i][j] = dist
distance_matrix[j][i] = dist
return distance_matrix
# Compute the distances to the depot
def compute_distance_depots(depot_x, depot_y, customers_x, customers_y):
nb_customers = len(customers_x)
distance_depots = [None] * nb_customers
for i in range(nb_customers):
dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
distance_depots[i] = dist
return distance_depots
def compute_dist(xi, xj, yi, yj):
return math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
if __name__ == '__main__':
if len(sys.argv) < 2:
print("Usage: python pdptw.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 pdptw.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.libpdptw instances\lc101.txt
- Compilation / Execution (Linux)
-
g++ pdptw.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o pdptw./pdptw instances/lc101.txt
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>
using namespace hexaly;
using namespace std;
class Pdptw {
public:
// Hexaly Optimizer
HexalyOptimizer optimizer;
// Number of customers
int nbCustomers;
// Capacity of the trucks
int truckCapacity;
// Latest allowed arrival to depot
int maxHorizon;
// Demand for each customer
vector<int> demandsData;
// Earliest arrival for each customer
vector<int> earliestStartData;
// Latest departure from each customer
vector<int> latestEndData;
// Service time for each customer
vector<int> serviceTimeData;
// Index for pickup for each node
vector<int> pickUpIndex;
// Index for delivery for each node
vector<int> deliveryIndex;
// Distance matrix between customers
vector<vector<double>> distMatrixData;
// Distances between customers and depot
vector<double> distDepotData;
// Number of trucks
int nbTrucks;
// Decision variables
vector<HxExpression> customersSequences;
// Are the trucks actually used
vector<HxExpression> trucksUsed;
// Cumulated lateness in the solution (must be 0 for the solution to be valid)
HxExpression totalLateness;
// Number of trucks used in the solution
HxExpression nbTrucksUsed;
// Distance traveled by all the trucks
HxExpression totalDistance;
Pdptw() {}
// Read instance data
void readInstance(const string& fileName) { readInputPdptw(fileName); }
void solve(int limit) {
// Declare the optimization model
HxModel model = optimizer.getModel();
// Sequence of customers visited by each truck
customersSequences.resize(nbTrucks);
for (int k = 0; k < nbTrucks; ++k) {
customersSequences[k] = model.listVar(nbCustomers);
}
// All customers must be visited by exactly one truck
model.constraint(model.partition(customersSequences.begin(), customersSequences.end()));
// Create Hexaly arrays to be able to access them with "at" operators
HxExpression demands = model.array(demandsData.begin(), demandsData.end());
HxExpression earliest = model.array(earliestStartData.begin(), earliestStartData.end());
HxExpression latest = model.array(latestEndData.begin(), latestEndData.end());
HxExpression serviceTime = model.array(serviceTimeData.begin(), serviceTimeData.end());
HxExpression distMatrix = model.array();
for (int n = 0; n < nbCustomers; ++n) {
distMatrix.addOperand(model.array(distMatrixData[n].begin(), distMatrixData[n].end()));
}
HxExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());
trucksUsed.resize(nbTrucks);
vector<HxExpression> distRoutes(nbTrucks), endTime(nbTrucks), homeLateness(nbTrucks), lateness(nbTrucks);
// Pickups and deliveries
HxExpression customersSequencesArray = model.array(customersSequences.begin(), customersSequences.end());
for (int i = 0; i < nbCustomers; ++i) {
if (pickUpIndex[i] == -1) {
HxExpression pickUpListIndex = model.find(customersSequencesArray, i);
HxExpression deliveryListIndex = model.find(customersSequencesArray, deliveryIndex[i]);
model.constraint(pickUpListIndex == deliveryListIndex);
HxExpression pickupList = model.at(customersSequencesArray, pickUpListIndex);
HxExpression deliveryList = model.at(customersSequencesArray, deliveryListIndex);
model.constraint(model.indexOf(pickupList, i) < model.indexOf(deliveryList, deliveryIndex[i]));
}
}
for (int k = 0; k < nbTrucks; ++k) {
HxExpression sequence = customersSequences[k];
HxExpression c = model.count(sequence);
// A truck is used if it visits at least one customer
trucksUsed[k] = c > 0;
// The quantity needed in each route must not exceed the truck capacity at any point in the sequence
HxExpression demandLambda = model.createLambdaFunction(
[&](HxExpression i, HxExpression prev) { return prev + demands[sequence[i]]; });
HxExpression routeQuantity = model.array(model.range(0, c), demandLambda, 0);
HxExpression quantityLambda =
model.createLambdaFunction([&](HxExpression i) { return routeQuantity[i] <= truckCapacity; });
model.constraint(model.and_(model.range(0, c), quantityLambda));
// Distance traveled by truck k
HxExpression distLambda = model.createLambdaFunction(
[&](HxExpression i) { return model.at(distMatrix, sequence[i - 1], sequence[i]); });
distRoutes[k] = model.sum(model.range(1, c), distLambda) +
model.iif(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);
// End of each visit
HxExpression endLambda = model.createLambdaFunction([&](HxExpression i, HxExpression prev) {
return model.max(earliest[sequence[i]],
model.iif(i == 0, distDepot[sequence[0]],
prev + model.at(distMatrix, sequence[i - 1], sequence[i]))) +
serviceTime[sequence[i]];
});
endTime[k] = model.array(model.range(0, c), endLambda, 0);
// Arriving home after max_horizon
homeLateness[k] =
model.iif(trucksUsed[k], model.max(0, endTime[k][c - 1] + distDepot[sequence[c - 1]] - maxHorizon), 0);
// Completing visit after latest_end
HxExpression lateLambda = model.createLambdaFunction(
[&](HxExpression i) { return model.max(0, endTime[k][i] - latest[sequence[i]]); });
lateness[k] = homeLateness[k] + model.sum(model.range(0, c), lateLambda);
}
// Total lateness
totalLateness = model.sum(lateness.begin(), lateness.end());
// Total number of trucks used
nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());
// Total distance traveled
totalDistance = model.round(100 * model.sum(distRoutes.begin(), distRoutes.end())) / 100;
// Objective: minimize the number of trucks used, then minimize the distance traveled
model.minimize(totalLateness);
model.minimize(nbTrucksUsed);
model.minimize(totalDistance);
model.close();
// Parameterize the optimizer
optimizer.getParam().setTimeLimit(limit);
optimizer.solve();
}
/* Write the solution in a file with the following format:
* - number of trucks used and total distance
* - for each truck the customers visited (omitting the start/end at the depot) */
void writeSolution(const string& fileName) {
ofstream outfile;
outfile.exceptions(ofstream::failbit | ofstream::badbit);
outfile.open(fileName.c_str());
outfile << nbTrucksUsed.getValue() << " " << totalDistance.getDoubleValue() << endl;
for (int k = 0; k < nbTrucks; ++k) {
if (trucksUsed[k].getValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +2 is to put it back in 2...nbCustomers+2
// as in the data files (1 being the depot)
HxCollection customersCollection = customersSequences[k].getCollectionValue();
for (int i = 0; i < customersCollection.count(); ++i) {
outfile << customersCollection[i] + 1 << " ";
}
outfile << endl;
}
}
private:
// The input files follow the "Li & Lim" format
void readInputPdptw(const string& fileName) {
ifstream infile(fileName.c_str());
if (!infile.is_open()) {
throw std::runtime_error("File cannot be opened.");
}
string str;
long dump;
int depotX, depotY;
vector<int> customersX;
vector<int> customersY;
infile >> nbTrucks;
infile >> truckCapacity;
infile >> dump;
infile >> dump;
infile >> depotX;
infile >> depotY;
infile >> dump;
infile >> dump;
infile >> maxHorizon;
infile >> dump;
infile >> dump;
infile >> dump;
while (infile >> dump) {
int cx, cy, demand, ready, due, service, pick, delivery;
infile >> cx;
infile >> cy;
infile >> demand;
infile >> ready;
infile >> due;
infile >> service;
infile >> pick;
infile >> delivery;
customersX.push_back(cx);
customersY.push_back(cy);
demandsData.push_back(demand);
earliestStartData.push_back(ready);
latestEndData.push_back(due + service); // in input files due date is meant as latest start time
serviceTimeData.push_back(service);
pickUpIndex.push_back(pick - 1);
deliveryIndex.push_back(delivery - 1);
}
nbCustomers = customersX.size();
computeDistanceMatrix(depotX, depotY, customersX, customersY);
infile.close();
}
// Compute the distance matrix
void computeDistanceMatrix(int depotX, int depotY, const vector<int>& customersX, const vector<int>& customersY) {
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) {
double distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
distMatrixData[i][j] = distance;
distMatrixData[j][i] = distance;
}
}
distDepotData.resize(nbCustomers);
for (int i = 0; i < nbCustomers; ++i) {
distDepotData[i] = computeDist(depotX, customersX[i], depotY, customersY[i]);
}
}
double computeDist(int xi, int xj, int yi, int yj) {
return sqrt(pow((double)xi - xj, 2) + pow((double)yi - yj, 2));
}
};
int main(int argc, char** argv) {
if (argc < 2) {
cerr << "Usage: pdptw 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 {
Pdptw 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 Pdptw.cs /reference:Hexaly.NET.dllPdptw instances\lc101.txt
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;
public class Pdptw : IDisposable
{
// Hexaly Optimizer
HexalyOptimizer optimizer;
// Number of customers
int nbCustomers;
// Capacity of the trucks
int truckCapacity;
// Latest allowed arrival to depot
int maxHorizon;
// Demand for each customer
List<int> demandsData;
// Earliest arrival for each customer
List<int> earliestStartData;
// Latest departure from each customer
List<int> latestEndData;
// Service time for each customer
List<int> serviceTimeData;
// Index for pick up for each node
List<int> pickUpIndex;
// Index for delivery for each node
List<int> deliveryIndex;
// Distance matrix between customers
double[][] distMatrixData;
// Distances between customers and depot
double[] distDepotData;
// Number of trucks
int nbTrucks;
// Decision variables
HxExpression[] customersSequences;
// Are the trucks actually used
HxExpression[] trucksUsed;
// Distance traveled by each truck
HxExpression[] distRoutes;
// End time array for each truck
HxExpression[] endTime;
// Home lateness for each truck
HxExpression[] homeLateness;
// Cumulated Lateness for each truck
HxExpression[] lateness;
// Cumulated lateness in the solution (must be 0 for the solution to be valid)
HxExpression totalLateness;
// Number of trucks used in the solution
HxExpression nbTrucksUsed;
// Distance traveled by all the trucks
HxExpression totalDistance;
public Pdptw()
{
optimizer = new HexalyOptimizer();
}
/* Read instance data */
void ReadInstance(string fileName)
{
ReadInputPdptw(fileName);
}
public void Dispose()
{
if (optimizer != null)
optimizer.Dispose();
}
void Solve(int limit)
{
// Declare the optimization model
HxModel model = optimizer.GetModel();
trucksUsed = new HxExpression[nbTrucks];
customersSequences = new HxExpression[nbTrucks];
distRoutes = new HxExpression[nbTrucks];
endTime = new HxExpression[nbTrucks];
homeLateness = new HxExpression[nbTrucks];
lateness = new HxExpression[nbTrucks];
// Sequence of customers visited by each truck
for (int k = 0; k < nbTrucks; ++k)
customersSequences[k] = model.List(nbCustomers);
// All customers must be visited by exactly one truck
model.Constraint(model.Partition(customersSequences));
// Create HexalyOptimizer arrays to be able to access them with "at" operators
HxExpression demands = model.Array(demandsData);
HxExpression earliest = model.Array(earliestStartData);
HxExpression latest = model.Array(latestEndData);
HxExpression serviceTime = model.Array(serviceTimeData);
HxExpression distDepot = model.Array(distDepotData);
HxExpression distMatrix = model.Array(distMatrixData);
HxExpression customersSequencesArray = model.Array(customersSequences);
for (int i = 0; i < nbCustomers; ++i)
{
if (pickUpIndex[i] == -1)
{
HxExpression pickUpListIndex = model.Find(customersSequencesArray, i);
HxExpression deliveryListIndex = model.Find(customersSequencesArray, deliveryIndex[i]);
model.Constraint(pickUpListIndex == deliveryListIndex);
HxExpression pickUpList = customersSequencesArray[pickUpListIndex];
HxExpression deliveryList = customersSequencesArray[deliveryListIndex];
model.Constraint(
model.IndexOf(pickUpList, i) < model.IndexOf(deliveryList, deliveryIndex[i])
);
}
}
for (int k = 0; k < nbTrucks; ++k)
{
HxExpression sequence = customersSequences[k];
HxExpression c = model.Count(sequence);
// A truck is used if it visits at least one customer
trucksUsed[k] = c > 0;
// The quantity needed in each route must not exceed the truck capacity at any point in the sequence
HxExpression demandLambda = model.LambdaFunction(
(i, prev) => prev + demands[sequence[i]]
);
HxExpression routeQuantity = model.Array(model.Range(0, c), demandLambda, 0);
HxExpression quantityLambda = model.LambdaFunction(
i => routeQuantity[i] <= truckCapacity
);
model.Constraint(model.And(model.Range(0, c), quantityLambda));
// Distance traveled by truck k
HxExpression distLambda = model.LambdaFunction(
i => distMatrix[sequence[i - 1], sequence[i]]
);
distRoutes[k] =
model.Sum(model.Range(1, c), distLambda)
+ model.If(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);
// End of each visit
HxExpression endLambda = model.LambdaFunction(
(i, prev) =>
model.Max(
earliest[sequence[i]],
model.If(
i == 0,
distDepot[sequence[0]],
prev + distMatrix[sequence[i - 1], sequence[i]]
)
) + serviceTime[sequence[i]]
);
endTime[k] = model.Array(model.Range(0, c), endLambda, 0);
// Arriving home after max_horizon
homeLateness[k] = model.If(
trucksUsed[k],
model.Max(0, endTime[k][c - 1] + distDepot[sequence[c - 1]] - maxHorizon),
0
);
// Completing visit after latest_end
HxExpression lateLambda = model.LambdaFunction(
i => model.Max(endTime[k][i] - latest[sequence[i]], 0)
);
lateness[k] = homeLateness[k] + model.Sum(model.Range(0, c), lateLambda);
}
totalLateness = model.Sum(lateness);
nbTrucksUsed = model.Sum(trucksUsed);
totalDistance = model.Round(100 * model.Sum(distRoutes)) / 100;
// Objective: minimize the number of trucks used, then minimize the distance traveled
model.Minimize(totalLateness);
model.Minimize(nbTrucksUsed);
model.Minimize(totalDistance);
model.Close();
// Parameterize the optimizer
optimizer.GetParam().SetTimeLimit(limit);
optimizer.Solve();
}
/* Write the solution in a file with the following format:
* - number of trucks used and total distance
* - for each truck the customers visited (omitting the start/end at the depot) */
void WriteSolution(string fileName)
{
using (StreamWriter output = new StreamWriter(fileName))
{
output.WriteLine(nbTrucksUsed.GetValue() + " " + totalDistance.GetDoubleValue());
for (int k = 0; k < nbTrucks; ++k)
{
if (trucksUsed[k].GetValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +2 is to put it back in 2...nbCustomers+2
// as in the data files (1 being the depot)
HxCollection customersCollection = customersSequences[k].GetCollectionValue();
for (int i = 0; i < customersCollection.Count(); ++i)
output.Write((customersCollection[i] + 1) + " ");
output.WriteLine();
}
}
}
public static void Main(string[] args)
{
if (args.Length < 1)
{
Console.WriteLine("Usage: Pdptw 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 (Pdptw model = new Pdptw())
{
model.ReadInstance(instanceFile);
model.Solve(int.Parse(strTimeLimit));
if (outputFile != null)
model.WriteSolution(outputFile);
}
}
private string[] SplitInput(StreamReader input)
{
string line = input.ReadLine();
if (line == null)
return new string[0];
return line.Split(new[] { '\t' }, StringSplitOptions.RemoveEmptyEntries);
}
// The input files follow the "Li & Lim" format
private void ReadInputPdptw(string fileName)
{
using (StreamReader input = new StreamReader(fileName))
{
string[] splitted;
splitted = SplitInput(input);
nbTrucks = int.Parse(splitted[0]);
truckCapacity = int.Parse(splitted[1]);
splitted = SplitInput(input);
int depotX = int.Parse(splitted[1]);
int depotY = int.Parse(splitted[2]);
maxHorizon = int.Parse(splitted[5]);
List<int> customersX = new List<int>();
List<int> customersY = new List<int>();
demandsData = new List<int>();
earliestStartData = new List<int>();
latestEndData = new List<int>();
serviceTimeData = new List<int>();
pickUpIndex = new List<int>();
deliveryIndex = new List<int>();
while (!input.EndOfStream)
{
splitted = SplitInput(input);
if (splitted.Length < 9)
break;
customersX.Add(int.Parse(splitted[1]));
customersY.Add(int.Parse(splitted[2]));
demandsData.Add(int.Parse(splitted[3]));
int ready = int.Parse(splitted[4]);
int due = int.Parse(splitted[5]);
int service = int.Parse(splitted[6]);
pickUpIndex.Add(int.Parse(splitted[7]) - 1);
deliveryIndex.Add(int.Parse(splitted[8]) - 1);
earliestStartData.Add(ready);
latestEndData.Add(due + service); // in input files due date is meant as latest start time
serviceTimeData.Add(service);
}
nbCustomers = customersX.Count;
ComputeDistanceMatrix(depotX, depotY, customersX, customersY);
}
}
// Compute the distance matrix
private void ComputeDistanceMatrix(
int depotX,
int depotY,
List<int> customersX,
List<int> customersY
)
{
distMatrixData = new double[nbCustomers][];
for (int i = 0; i < nbCustomers; ++i)
distMatrixData[i] = new double[nbCustomers];
for (int i = 0; i < nbCustomers; ++i)
{
distMatrixData[i][i] = 0;
for (int j = i + 1; j < nbCustomers; ++j)
{
double dist = ComputeDist(
customersX[i],
customersX[j],
customersY[i],
customersY[j]
);
distMatrixData[i][j] = dist;
distMatrixData[j][i] = dist;
}
}
distDepotData = new double[nbCustomers];
for (int i = 0; i < nbCustomers; ++i)
distDepotData[i] = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
}
private double ComputeDist(int xi, int xj, int yi, int yj)
{
return Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
}
}
- Compilation / Execution (Windows)
-
javac Pdptw.java -cp %HX_HOME%\bin\hexaly.jarjava -cp %HX_HOME%\bin\hexaly.jar;. Pdptw instances\lc101.txt
- Compilation / Execution (Linux)
-
javac Pdptw.java -cp /opt/hexaly_13_0/bin/hexaly.jarjava -cp /opt/hexaly_13_0/bin/hexaly.jar:. Pdptw instances/lc101.txt
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;
public class Pdptw {
// Hexaly Optimizer
private final HexalyOptimizer optimizer;
// Number of customers
int nbCustomers;
// Capacity of the trucks
private int truckCapacity;
// Latest allowed arrival to depot
int maxHorizon;
// Demand for each customer
List<Integer> demandsData;
// Earliest arrival for each customer
List<Integer> earliestStartData;
// Latest departure from each customer
List<Integer> latestEndData;
// Service time for each customer
List<Integer> serviceTimeData;
// Index for pick up for each node
List<Integer> pickUpIndex;
// Index for delivery for each node
List<Integer> deliveryIndex;
// Distance matrix between customers
private double[][] distMatrixData;
// Distances between customers and depot
private double[] distDepotData;
// Number of trucks
private int nbTrucks;
// Decision variables
private HxExpression[] customersSequences;
// Are the trucks actually used
private HxExpression[] trucksUsed;
// Distance traveled by each truck
private HxExpression[] distRoutes;
// End time array for each truck
private HxExpression[] endTime;
// Home lateness for each truck
private HxExpression[] homeLateness;
// Cumulated Lateness for each truck
private HxExpression[] lateness;
// Cumulated lateness in the solution (must be 0 for the solution to be valid)
private HxExpression totalLateness;
// Number of trucks used in the solution
private HxExpression nbTrucksUsed;
// Distance traveled by all the trucks
private HxExpression totalDistance;
private Pdptw() {
optimizer = new HexalyOptimizer();
}
// Read instance data
private void readInstance(String fileName) throws IOException {
readInputPdptw(fileName);
}
private void solve(int limit) {
// Declare the optimization model
HxModel m = optimizer.getModel();
trucksUsed = new HxExpression[nbTrucks];
customersSequences = new HxExpression[nbTrucks];
distRoutes = new HxExpression[nbTrucks];
endTime = new HxExpression[nbTrucks];
homeLateness = new HxExpression[nbTrucks];
lateness = new HxExpression[nbTrucks];
// Sequence of customers visited by each truck
for (int k = 0; k < nbTrucks; ++k)
customersSequences[k] = m.listVar(nbCustomers);
// All customers must be visited by exactly one truck
m.constraint(m.partition(customersSequences));
// Create HexalyOptimizer arrays to be able to access them with "at" operators
HxExpression demands = m.array(demandsData);
HxExpression earliest = m.array(earliestStartData);
HxExpression latest = m.array(latestEndData);
HxExpression serviceTime = m.array(serviceTimeData);
HxExpression distDepot = m.array(distDepotData);
HxExpression distMatrix = m.array(distMatrixData);
// Pickups and deliveries
HxExpression customersSequencesArray = m.array(customersSequences);
for (int i = 0; i < nbCustomers; ++i) {
if (pickUpIndex.get(i) == -1) {
HxExpression pickUpListIndex = m.find(customersSequencesArray, i);
HxExpression deliveryListIndex = m.find(customersSequencesArray, deliveryIndex.get(i));
m.constraint(m.eq(pickUpListIndex, deliveryListIndex));
HxExpression pickUpList = m.at(customersSequencesArray, pickUpListIndex);
HxExpression deliveryList = m.at(customersSequencesArray, deliveryListIndex);
m.constraint(m.leq(m.indexOf(pickUpList, i), m.indexOf(deliveryList, deliveryIndex.get(i))));
}
}
for (int k = 0; k < nbTrucks; ++k) {
HxExpression sequence = customersSequences[k];
HxExpression c = m.count(sequence);
// A truck is used if it visits at least one customer
trucksUsed[k] = m.gt(c, 0);
// The quantity needed in each route must not exceed the truck capacity at any point in the sequence
HxExpression demandLambda = m.lambdaFunction((i, prev) -> m.sum(prev, m.at(demands, m.at(sequence, i))));
HxExpression routeQuantity = m.array(m.range(0, c), demandLambda, 0);
HxExpression quantityLambda = m.lambdaFunction(i -> m.leq(m.at(routeQuantity, i), truckCapacity));
m.constraint(m.and(m.range(0, c), quantityLambda));
// Distance traveled by truck k
HxExpression distLambda = m
.lambdaFunction(i -> m.at(distMatrix, m.at(sequence, m.sub(i, 1)), m.at(sequence, i)));
distRoutes[k] = m.sum(m.sum(m.range(1, c), distLambda), m.iif(m.gt(c, 0),
m.sum(m.at(distDepot, m.at(sequence, 0)), m.at(distDepot, m.at(sequence, m.sub(c, 1)))), 0));
// End of each visit
HxExpression endLambda = m.lambdaFunction((i, prev) -> m.sum(
m.max(m.at(earliest, m.at(sequence, i)),
m.sum(m.iif(m.eq(i, 0), m.at(distDepot, m.at(sequence, 0)),
m.sum(prev, m.at(distMatrix, m.at(sequence, m.sub(i, 1)), m.at(sequence, i)))))),
m.at(serviceTime, m.at(sequence, i))));
endTime[k] = m.array(m.range(0, c), endLambda, 0);
HxExpression theEnd = endTime[k];
// Arriving home after max_horizon
homeLateness[k] = m.iif(trucksUsed[k],
m.max(0,
m.sum(m.at(theEnd, m.sub(c, 1)), m.sub(m.at(distDepot, m.at(sequence, m.sub(c, 1))), maxHorizon))),
0);
// Completing visit after latest_end
HxExpression lateLambda = m
.lambdaFunction(i -> m.max(m.sub(m.at(theEnd, i), m.at(latest, m.at(sequence, i))), 0));
lateness[k] = m.sum(homeLateness[k], m.sum(m.range(0, c), lateLambda));
}
totalLateness = m.sum(lateness);
nbTrucksUsed = m.sum(trucksUsed);
totalDistance = m.div(m.round(m.prod(100, m.sum(distRoutes))), 100);
// Objective: minimize the number of trucks used, then minimize the distance traveled
m.minimize(totalLateness);
m.minimize(nbTrucksUsed);
m.minimize(totalDistance);
m.close();
// Parameterize the optimizer
optimizer.getParam().setTimeLimit(limit);
optimizer.solve();
}
/*
* Write the solution in a file with the following format:
* - number of trucks used and total distance
* - for each truck the customers visited (omitting the start/end at the depot)
*/
private void writeSolution(String fileName) throws IOException {
try (PrintWriter output = new PrintWriter(fileName)) {
output.println(nbTrucksUsed.getValue() + " " + totalDistance.getDoubleValue());
for (int k = 0; k < nbTrucks; ++k) {
if (trucksUsed[k].getValue() != 1)
continue;
// Values in sequence are in 0...nbCustomers. +2 is to put it back in 2...nbCustomers+2
// as in the data files (1 being the depot)
HxCollection customersCollection = customersSequences[k].getCollectionValue();
for (int i = 0; i < customersCollection.count(); ++i) {
output.print((customersCollection.get(i) + 1) + " ");
}
output.println();
}
}
}
// The input files follow the "Li & Lim" format
private void readInputPdptw(String fileName) throws IOException {
try (Scanner input = new Scanner(new File(fileName))) {
nbTrucks = input.nextInt();
truckCapacity = input.nextInt();
input.nextInt();
input.nextInt();
int depotX = input.nextInt();
int depotY = input.nextInt();
input.nextInt();
input.nextInt();
maxHorizon = input.nextInt();
input.nextInt();
input.nextInt();
input.nextInt();
List<Integer> customersX = new ArrayList<Integer>();
List<Integer> customersY = new ArrayList<Integer>();
demandsData = new ArrayList<Integer>();
earliestStartData = new ArrayList<Integer>();
latestEndData = new ArrayList<Integer>();
serviceTimeData = new ArrayList<Integer>();
pickUpIndex = new ArrayList<Integer>();
deliveryIndex = new ArrayList<Integer>();
while (input.hasNextInt()) {
input.nextInt();
int cx = input.nextInt();
int cy = input.nextInt();
int demand = input.nextInt();
int ready = input.nextInt();
int due = input.nextInt();
int service = input.nextInt();
int pick = input.nextInt();
int delivery = input.nextInt();
customersX.add(cx);
customersY.add(cy);
demandsData.add(demand);
earliestStartData.add(ready);
latestEndData.add(due + service); // in input files due date is meant as latest start time
serviceTimeData.add(service);
pickUpIndex.add(pick - 1);
deliveryIndex.add(delivery - 1);
}
nbCustomers = customersX.size();
computeDistanceMatrix(depotX, depotY, customersX, customersY);
}
}
// Compute the distance matrix
private void computeDistanceMatrix(int depotX, int depotY, List<Integer> customersX, List<Integer> customersY) {
distMatrixData = new double[nbCustomers][nbCustomers];
for (int i = 0; i < nbCustomers; ++i) {
distMatrixData[i][i] = 0;
for (int j = i + 1; j < nbCustomers; ++j) {
double dist = computeDist(customersX.get(i), customersX.get(j), customersY.get(i), customersY.get(j));
distMatrixData[i][j] = dist;
distMatrixData[j][i] = dist;
}
}
distDepotData = new double[nbCustomers];
for (int i = 0; i < nbCustomers; ++i) {
distDepotData[i] = computeDist(depotX, customersX.get(i), depotY, customersY.get(i));
}
}
private double computeDist(int xi, int xj, int yi, int yj) {
return Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
}
public static void main(String[] args) {
if (args.length < 1) {
System.err.println("Usage: java Pdptw inputFile [outputFile] [timeLimit]");
System.exit(1);
}
try {
String instanceFile = args[0];
String outputFile = args.length > 1 ? args[1] : null;
String strTimeLimit = args.length > 2 ? args[2] : "20";
Pdptw model = new Pdptw();
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 reaches an average gap of 1.1% on the Pickup and Delivery Problem with Time Windows (PDPTW) in 1 minute of running time, on the instances from the Li & Lim benchmark, with 100 to 1,000 customers. Our Pickup and Delivery Problem with Time Windows (PDPTW) benchmark page illustrates how Hexaly Optimizer outperforms traditional general-purpose optimization solvers like OR-Tools on this fundamental but challenging problem.