key: cord-0861815-ikmfsw2k
authors: Dzalbs, Ivars; Kalganova, Tatiana
title: Accelerating supply chains with Ant Colony Optimization across a range of hardware solutions
date: 2020-06-29
journal: Comput Ind Eng
DOI: 10.1016/j.cie.2020.106610
sha: 4d74e5ae710c0fc3d62a8c8aec4a79c6961867e7
doc_id: 861815
cord_uid: ikmfsw2k

Ant Colony algorithm has been applied to various optimisation problems, however, most of the previous work on scaling and parallelism focuses on Travelling Salesman Problems (TSPs). Although useful for benchmarks and new idea comparison, the algorithmic dynamics do not always transfer to complex real-life problems, where additional meta-data is required during solution construction. This paper explores how the benchmark performance differs from real-world problems in the context of Ant Colony Optimization (ACO) and demonstrate that in order to generalise the findings, the algorithms have to be tested on both standard benchmarks and real-world applications. ACO and its scaling dynamics with two parallel ACO architectures – Independent Ant Colonies (IAC) and Parallel Ants (PA). Results showed that PA was able to reach a higher solution quality in fewer iterations as the number of parallel instances increased. Furthermore, speed performance was measured across three different hardware solutions – 16 core CPU, 68 core Xeon Phi and up to 4 Geforce GPUs. State of the art, ACO vectorisation techniques such as SS-Roulette were implemented using C++ and CUDA. Although excellent for routing simple TSPs, it was concluded that for complex real-world supply chain routing GPUs are not suitable due to meta-data access footprint required. Thus, our work demonstrates that the standard benchmarks are not suitable for generalised conclusions.

Supply chain optimisation has become an integral part of any global company with a complex manufacturing and distribution network. For many companies, inefficient distribution plan can make a significant difference to the bottom line. Modelling a complete distribution network from the initial materials to the delivery to the customer is very computationally intensive. With increasing supply chain modelling complexity in ever-changing global geo-political environment, fast adaptability is an edge. A company can model the impact of currency exchange rate changes, import tax policy reforms, oil price fluctuations and political events such as Brexit, Covid-19 before they happen. Such modelling requires fast optimisation algorithms. Mixed Integer Linear Programming (MILP) tools such as Cplex are commonly used to optimise various supply chain networks [1] . Although MILP tools can obtain an optimum solution for a large variety of linear models, not all real-world supply chain models are linear. Furthermore, MILP is computationally expensive and on large instances can fail to produce an optimal solution. For that reason, many alternative algorithmic approaches (heuristics, meta-heuristics, fuzzy methods) have been explored to solve large-complex SC models [1] . One of these algorithms is the Ant Colony Optimization (ACO), which can be well mapped to real-world problems such as routing [2] and scheduling [3] . Supply Chain Optimization Problem (SCOP) includes both, finding the best route to ship a specific order and finding the most optimal time to ship it, such that it reaches expected customer satisfaction while minimising the total cost occurred. Although other metaheuristics algorithms exist in the literature for solving SCOPs, such as Genetic Algorithm (GA) [4] [5] and Simulated Annealing (SA) [6] [7] , ACO was chosen due to the long history of the algorithm applied to various vehicle routing [8] [9] and supply chain [10] [11] problems with great solution quality and speed. Also, a recent study in [12] concluded that compared to GA and SA, the ACO performs the best for routing problems such as the Travelling Salesman Problem (TSP). Ant colony algorithms try to mimic the observed behaviour of ants inside colonies to solve a large range of optimisation problems. Since the introduction by Marco Dorigo in 1992, many variations and hybrid approaches of Ant Colony algorithms have been explored [13] [14] . Most ant colony algorithms consist of two distinct stages -solution construction and pheromone feedback to other ants. Typically, an artificial ant builds its solution from the pheromone left from previous ants, therefore allowing communication over many iterations via a pheromone matrix and converges to a better solution. The process of solution creation and pheromone update is repeated over many iterations until the termination criterion is reached; this can be either total number of iterations, total computation time or dynamic termination. Researchers in [15] compared an industrial optimisation-based tool -IBM ILOG Cplex with their proposed ACO algorithm. It was concluded that the proposed algorithm covered 94% of optimal solutions on small problems and 88% for large-size problems while consuming significantly less computation time. Similarly, [16] compared ACO and Cplex performance on multi-product and multi-period Inventory Routing Problem. On small instances, ACO reached 95% of the optimal solution while on large instances performed better than time-constrained Cplex solver. Furthermore, ACO implementations of Closed-Loop Supply Chain (CLSC) have been proposed; CLSC contains two parts of the supply chain -forward supply and reverse/return. [17] solved CLSC models, where the ACO implementation outperformed commercial MILP (Cplex) on nonlinear instances and obtained 98% optimal solution with 40% less computation time on linear instances.

Academic literature suggests that Graphical Processing Units (GPUs) are very suitable for solving benchmark routing problems such as Travelling Salesman Problem (TSP), with speedups of up to 60x [18] and even 172x [19] when compared to the sequential CPU implementation. This paper aims to explore if the same ACO architectures that are so well suited for TSP can be applied for a real-world supply chain optimisation problem. Furthermore, investigate what hardware architectures are the best suited for the supply chain problem solved. The paper is structured as follows: Section 2 explores the current state of the art parallel implementations of ACO across CPU, GPU and Xeon Phi; Section 3 introduces the hardware and software solutions used; Section 4 described the real-world problem being solved; Section 5 details the parallel ACO implementations and Section 6 compares the results. Finally, Section 7 concludes the paper.

Since the introduction of ACO in 1992, numerous ACO algorithms have been applied to many different problems, and many different parallel architectures have been explored previously. [20] specifies 5 of such architectures:

 Parallel Independent Ant Colonies -each ant colony develop their solutions in parallel without any communication in-between;  Parallel Interacting Ant Colonies -each colony creates a solution in parallel and some information is shared between the colonies;  Parallel Ants -each ant builds solution independently, then all the resulting pheromones are shared for the next iteration;  Parallel Evaluation of Solution Elements -for problems where fitness function calculations take considerably more time than the solution creation;  Parallel Combination of Ants and Evaluation of Solution Elements -a combination of any of the above.

Researchers have tried to exploit the parallelism offered from recent multi-core CPUs [21] , along with clusters of CPUs ( [22] [23] ) and most recently GPUs [24] and Intel's many-core architectures such as Xeon Phi [25] . Breakdown of the strategies and problems solved are shown in Table 1 . Table 2 shows the most common problems solved by ACO and their corresponding associated constraints / meta-data required during solution creation. 

Parallel ACO CPU architectures have been applied to various tasks -for example, [26] applied ACO for supply chain scheduling problem in mining domain. Authors managed to reduce the execution time from one hour (serial) to around 7 minutes. Both [40] and [41] used ACO for image edge detection with varying results, [40] achieved a speedup of 3-5 times while [41] managed to reduce sequential runtime by 30%. Most commonly, ACO has been applied to the Travelling Salesman Problem (TSP) benchmarks. For instance, [28] proposed ACO approach with randomly synchronised ants, the strategy showed a faster convergence compared to other TSP approaches. Moreover, authors in [29] proposed a new multi-core Single Instruction Multiple Data (SIMD) model for solving TSPs. Similarly, both [42] and [43] tries to solve large instances of TSP (up to 200k and 20k cities, respectively) where the architectures are limited to the size of the pheromone matrix.

[43] discusses such limitations and proposes a new pheromone sharing for local searcheffective heuristics ACO (ESACO), which was able to compute TSP instances of 20k. In contrast, authors in [42] eliminate the need for pheromone matrix and store only the best solutions similar to the Population ACO. Furthermore, researchers implement a Partial Ant, also known as the cunning ant, where ant takes an existing partial solution and builds on top of it. Speedups of as much as 1200x are achieved compared to sequential Population ACO. Generally, CPU parallel architecture implementations come down to three programming approaches -Message Passing Interface (MPI) parallelism, OpenMP parallelism [44] and data parallelism with the vectorisation of SIMD. For instance, [45] explored both masterslave and coarse-grained strategies for ACO parallelisation using MPI. It was concluded that fine-grained master-slave policy performed the best. [46] used MPI with ACO to accelerate Maximum Weight Clique Problem (MWCP). The proposed algorithm was comparable to the ones in literature and outperformed Cplex solver in both -time and performance. Moreover, authors in [36] implemented parallel ACO for solving Flow shop scheduling problem with restrictions using MPI. Compared to the sequential version of the algorithm, 93 node cluster achieved a speedup of 16x. [47] compared ACO parallel implementation on MPI and OpenMP on small vector estimation problem. It was found that maximum speedup of OpenMP was 24x while MPI -16x. Furthermore, [29] explored the multi-core SIMD CPU with OpenCL and compared it to the performance of the GPU. It was found optimised parallel CPU-SIMD version can achieve similar solution quality and computation time than the state of art GPU implementation solving TSP.

Intel's Xeon Phi Many Integrated Core (MIC) architecture offers many cores on the CPU (60-72 cores per node) while offering lower clock frequency. Few researchers have had the opportunity to research ACO on the Xeon Phi architecture. For instance, [37] showed how utilising L1 and L2 cache on Xeon Phi coprocessor allowed a speedup of 42x solving TSP compared to sequential execution. Due to the nature of SIMD features such as AVX-512 on Xeon Phi, researchers in both [38] and [39] proposed a vectorisation model for roulette wheel selection in TSP. In the case of [39] , a 16.6x speedup was achieved compared to sequential execution. To the best of the authors' knowledge, Xeon Phi and ACO parallelism have not been explored to any other problem except TSP.

General Purpose GPU (GPGPU) programming is a growing field in computer science and machine learning. Many researchers have tried exploiting latest GPU architectures to speed optimise the convergence of ACO. ACO GPU implementation expands to many fields, such as edge detection ([35] [48] ), protein folding [30] , solving Multidimensional Knapsack Problems (MKPs) [31] and Vertex colouring problems [49] . Although task parallelism has potential for a speedup, [62] showed how data parallelism (vectorisation) on GPU could achieve better performance by proposed Independent Roulette wheel (I-Roulette). Authors then expanded the I-Roulette implementation in [33] , where SS-Roulette wheel was introduced. SS-Roulette stands for Scan and Stencil Roulette wheel. It mimics a sequential roulette wheel while allowing higher throughput due to parallelism. First, the Tabu list is multiplied by the probabilities and the results stored in a choice vector (scan). Then a stencil pattern is applied to the choice vector based on a random number to select an individual (stencil). Further, [19] implements a G-Roulette -a grouped roulette wheel selection based on I-Roulette, where cities in TSP selection are grouped in CUDA warps 1 . An impressive speedup of 172x was achieved compared to the sequential counterpart.

Fairly comparing parallel performances of different hardware architectures is by no means trivial. Most research compares a sequential CPU ACO implementation to one of the parallel GPUs, which is hardly fair [63] . In addition, unoptimised sequential code is compared to highly optimised GPU code. Such comparisons result in misleading and inflated speedups [24] . Furthermore, [32] argues that often the parameter settings chosen for the sequential implementation is biased in favour of GPU. [24] proposes criteria to calculate the real-world efficiency of two different hardware architectures by comparing the theoretical peak performances of GPU and CPU. While the proposed method is more appropriate, it still does not account for real-life scenarios where memory latency/speed, cache size, compilers and operating systems all play a role of the final execution time.

Therefore, two different systems with similar theoretical floating-point operations per second running the same executable can have significantly different execution times. Furthermore, in some instances, only execution time or solution quality is compared, rarely both are taken into consideration when analysing results.

This section briefly covers the tools and hardware-specific languages used in the implementation.

OpenMP 2 is a set of directives to a compiler that allows a programmer to create parallel tasks as well as vectorisation (Single Instruction Multiple Data -SIMD) to speed up execution of a program. A program containing parallel OpenMP directives starts as a single thread. Once directive such as #pragma omp parallel is reached, the main thread will create a thread pool and all methods within the #pragma region will be executed in parallel by each thread in the thread group. Moreover, once the thread reaches the end of the region, it will wait for all other threads to finish before dissolving the thread group and only the main thread will continue. Furthermore, OpenMP also supports nesting, meaning a thread in a thread-group can create its own individual thread-group and become the master thread for the newly created thread-group. However, thread-group creation and elimination can have significant overhead and therefore, thread-group re-use is highly recommended [64] . This paper utilises both omp parallel and omp simd directives.

Compute Unified Device Architecture (CUDA) is a General-purpose computing model on GPU developed by Nvidia in 2006. Since then, this proprietary framework has been utilised in the high-performance computing space via multiple Artificial Intelligence (AI) and Machine Learning (ML) interfaces and libraries/APIs. CUDA allows writing C programs that take advantage of any recent Nvidia GPU found in laptops, workstations and data centres. Each GPU contains multiple Streaming Multiprocessors (SM) that are designed to execute hundreds of threads concurrently. To achieve that, CUDA implements SIMT (Single Instruction Multiple-Threads) architecture, where instructions are pipelined for instructionlevel parallelism. Threads are grouped in sets of 32 -called warps. Each warp executes one instruction at a time on each thread. Furthermore, CUDA threads can access multiple memory spaces -global memory (large size, slower), texture memory (read only), shared memory (shared across threads in the same SM, lower latency) and local memory (limited set of registers within each thread, fastest) 3 . A batch of threads is grouped into a thread-block. Multiple thread-blocks create a grid of thread blocks. The programmer specifies the grid dimensionality at kernel launch time, by providing the number of thread-blocks and the number of threads per thread-block. Kernel launch fails if the program exceeds the hardware resource boundaries.

Knights Landing is a product code name for Intel's second-generation Intel Xeon Phi processors. First-generation of Xeon Phi, named Knights Corner, was a PCI-e coprocessor card based on many Intel Atom processor cores and support for Vector Processing Units (VPUs). The main advancement over Knights Corner was the standalone processor that can boot stock operating systems, along with improved power efficiency and vector performance. Furthermore, it also introduced a new high bandwidth Multi-Channel DRAM (MCDRAM) memory. Xeon phi support for standard x86 and x86-64 instructions, allows majority CPU compiled binaries to run without any modification. Moreover, support for 512-bit Advanced Vector Extensions (AVX-512) allows high throughput vector manipulations. Figure 1 . Knights Landing tile with a larger processor die [65] The Knights Landing cores are divided into tiles (typically between 32 and 36 tiles in total). Each tile contains two processor cores and each core is connected to two vector processing units (VPUs). Utilising AVX-512 and two VPUs, each core can deliver 32 dual-precision (DP) or 64 single-precision (SP) operations each cycle [65] . Furthermore, each core supports up to four threads of execution -hyper threads where instructions are pipelined.

Another introduction with the Knights Landing is the cluster modes and MCDRAM/DRAM management. The processor offers three primary cluster modes 4 -All to all mode, Quadrant mode and Sub-Numa Cluster (SNC) mode and three memory modes -cache mode, flat mode and hybrid mode. For a detailed description of the Knights Landing Xeon Phi architecture refer to [65] .

A real-world dataset of an outbound logistics network is provided by a global microchip producer. The company provided demand data for 9216 orders that need to be routed via their outbound supply chain network of 15 warehouses, 11 origin ports and one destination port (see Figure 2 ). Warehouses are limited to a specific set of products that they stock, furthermore, some warehouses are dedicated for supporting only a particular set of customers. Moreover, warehouses are limited by the number of orders that can be processed in a single day. A customer making an order decides what sort of service level they require -DTD (Door to Door), DTP (Door to Port) or CRF (Customer Referred Freight). In the case of CRF, the customer arranges the freight and company only incur the warehouse cost. In most instances, an order can be shipped via one of 9 couriers offering different rates for different weight bands and service levels. Although most of the shipments are made via air transport, some orders are shipped via ground -by trucks. The majority of couriers offer discounted rates as the total shipping weight increases based on different weight bands. However, a minimum charge for shipment still applies. Furthermore, faster shipping tends to be more expensive, but offer better customer satisfaction. Customer service level is out of the scope of this research. can be supplied by either origin ports or . In contrast, warehouse can only be

supplied via origin port and warehouse can be only supplied by origin port . In

VmiCustomers contains all edge cases, where the warehouse is only allowed to support specific customer, while any other non-listed warehouse can supply any customer. Moreover, the WhCapacities lists warehouse capacities measured in the number of orders per day and the WhCosts specifies the cost associated in storing the products in a given warehouse measured in dollars per unit.

The main goal of optimisation is to find a set of warehouses, shipping lanes and couriers to use for the most cost-effective supply chain. Therefore the fitness function is derived from two incurred costs -warehouse cost and transportation cost in equation (1). The totalling cost is then calculated across all orders in the dataset.

is warehouse cost for order at warehouse and is transportation cost for order between warehouse port and customer port , the total number of orders .

(2) Where warehouse cost for order at warehouse is calculated in (2), by the number of units in order multiplied by the warehouse storage rate (WhCosts table) .

Furthermore, transportation cost for a given order k and chosen line between origin port p and destination port j is calculated by the algorithm in Figure 3 : where is the service level for order , -origin port, -destination port, -courier, -service level, -delivery time, -transportation mode. Furthermore, is the minimum charge for given line , is the weight in kilograms for order . is the freight rate (dollars per kilogram) for given weight gap based on the total weight for the line (FreightRates table) .

The algorithm first checks what kind of service level the order requires, if the service level is equal to CRF (Customer Referred Freight) -transportation cost is 0. Furthermore, if order transportation mode is equal to GROUND (order transported via truck), order transportation cost is proportional to the weight consumed by the order ( ) in respect of the total weight for given line and the rate charged by a courier for full track . In all other cases, the transportation cost is calculated based on order weight and the freight rate . The freight rate is determined based on total weight on any given line and the corresponding weight band in the freight rate table. Furthermore, a minimum charge is applied in cases where the air transportation cost is less than the minimum charge.

The problem being solved complies with the following constraints:

Where = 1 if order was shipped from warehouse and 0 otherwise. is the order limit per day for warehouse (WhCapacities table) .

Where is the weight in kilograms for order shipped from warehouse port to customer port via courier using service level , delivery time and transportation mode . is the upper weight gap limit for line (FreightRates table) .

Where product for order belongs to supported products at warehouse (ProductsPerPlant table) . Warehouses can only support given customer in the VmiCustomers table, while all other warehouses that are not in the table can supply any customer. Moreover, the warehouse can only ship orders via supported origin port, defined in PlantPorts table.

To solve the transportation network optimisation problem, we are using an Ant Colony System algorithm first proposed by [67] . Because ACO is an iterative algorithm, it does require sequential execution. Therefore, the most naïve approach for parallel ACO is running multiple Independent Ant Colonies (IAC) with a unique seed for the pseudorandom number generator for each colony (high-level pseudocode in Figure 4 ). Due to the stochastic nature of solution creation, it is, therefore, more probabilistic to reach a better solution than a single colony. This approach has the advantage of low overhead as it requires no synchronisation between the parallel instances during the search. At the very end of the search, the best solution of all parallel colonies is chosen as the final solution. The main disadvantage of IAC is that if one of the colonies finds a better solution, there is no way to improve all the other colony's fitness values.

1. for all parallel instances m parallel do 2.

for all iterations iter do 3.

for all local ants a do 4.

local pheromone = global pheromone 5.

construct solution 6.

local pheromone update 7.

end for 8.

update global pheromone update based on the best solution 9.

end for 10. end for 11. find the best solution across parallel instances Alternatively, the ACO search algorithm could also be letting the artificial ant colonies synchronise after every iteration. Therefore, all parallel instances are aware of the best solution and can share pheromones accordingly. High-level pseudocode of such Parallel Ant (PA) implementation is shown in Figure 5 . The main advantage of this architecture is that it allows efficient pheromone sharing, therefore converging faster. However, there is a high risk of getting stuck into local optima as all ants start iteration with the same pheromone matrix. Furthermore, synchronisation of all parallel instances after every iteration is costly.

1. for all iterations iter do 2.

for all parallel instances m parallel do 3.

for all local ants a do 4.

local pheromone = global pheromone 5.

construct solution 6.

local pheromone update 7.

end for 8. end for 9.

find the best solution across parallel instances 10. update global pheromone update based on the best solution 11. end for Both IAC and PA implementations are exploiting task parallelism -each parallel instance (thread) gets a set of tasks to complete. An alternative approach would be to look at data parallelism and vectorisation. In such a strategy, each thread processes a specific section of the data and cooperatively complete the given task. Due to the highly sequential parts of ACO, it would not be practical to only use vectorisation alone. A more desirable path would be to implement vectorisation in conjugate to the task parallelism. In case of CPU, task parallelism can be done by the threads, while vectorisation is done by Vector Processing Units (VPUs) based on Advanced Vector Extensions 2 (AVX2) or AVX512. Moreover, in the case of GPU and CUDA -task parallelism would be done at a threadblock level while data parallelism would exploit WARP structures. Parallel Ants with Vectorisation (PAwV) expands on the Parallel Ants architecture by introducing dataparallelism of solution creation and an alternative roulette wheel implementation -SS-Roulette, first proposed in [68] . Local search in Figure 6 expands on the implementation in Figure 5 (lines 3-7) . First, the choiceMatrix is calculated by multiplying the probability of the route to be chosen with the tabuList -a list of still available routes (where 0 represents not available and 1 -route still can be selected). A random number between 0 and 1 is generated to determine if a given route will be chosen based on exploitation or exploration. In the case of exploitation, the choiceMatrix is reduced to obtain the maximum and the corresponding route index. Furthermore, in the case of exploration, the route is chosen based on the SS-Roulette wheel described by [68] . The main advantage of IAC is that it requires to synchronise between threads only at the start of the search and at the very end of the search, therefore keeping synchronisation overhead low. However, as there is no pheromone sharing, new better solutions cannot be shared across the parallel instances. In contrast, both PA and PAwV offers sharing of the best performing ants' pheromone before the next iteration begins. The potential drawback is that search might get stuck in local optimum as all parallel instances share the same pheromone starting point. Furthermore, pheromone sharing and therefore, synchronisation between threads is costly overhead, especially if performed after each iteration. The PAwV architecture exploits the use of SIMD instructions for further data parallelism inside the Ant's solution construction. Table 3 summarises these architectural features. 

A sequential implementation of ACO described in [67] is adapted from [69] by altering the heuristic information calculation for a given route -defined as a proportion of order's weight and the maximum weight gap (see Equation (2)). Furthermore, the ACO set of parameters were obtained from both work in [69] and empirical experimentation. Table 4 summarises these algorithm hyperparameters. Moreover, we then implement three different Parallel ACO architectures -Independent Ant Colonies (IAC), Parallel Ants (PA) and Parallel Ants with Vectorisation (PAwV) in C++ and CUDA C. Experiments were conducted on three different hardware configurations -CPU, GPU and Xeon Phi. Hardware architecture C shares the same host CPU as Hardware A.

It is crucial to consider both elapsed time and solution quality when referring to speed optimisation of optimisation algorithms. One could get superior convergence within iteration but, take twice as long to compute. Similarly, one could claim that the algorithm is much faster at completing a defined number of iterations but sacrifice solution quality. Furthermore, there is little point comparing sequential execution of one hardware platform to parallel implementation of another. A comparison should take into consideration all platform strengths and weaknesses and set up the most suitable configuration for a given platform. To obtain a baseline fitness convergence rate at a various number of parallel instances, we create Iterations vs Parallel Instances matrix for all architectures. An example of such matrix for Parallel Ants is shown in We then compute the number of iterations required to reach a specific solution quality for different ACO architectures in Table 6 , expressed as proximity to the best-known optimal solution. For the particular problem and dataset, the best solution is the total cost of 2,701,367.58. There are six checkpoints of solution quality ranging from 99% to 99.9%. Although at first 1% gain might not seem significant, one must remember that global supply chain costs are measured in hundreds of millions, and even 1% savings do affect the bottom line. Empty fields (-) represent instances where the ACO was not able to converge to given solution quality. On all experiments, IAC was able to obtain solution quality only below 99.6%. In contrast, PA and PA with 5 ant local search were able to achieve above 99.9% solution quality with 512 and 1024 parallel instances. Furthermore, IAC did not see any significant benefit of adding more parallel instances for 99% and 99.25% checkpoints. Table 5 . The matrix is derived by averaging the resulting fitness obtained from 10 independent simulations with a unique seed value for each given Parallel Instances configuration. All configurations are run for x number of iterations, where x is based on the total number of solutions explored and is a function of the number of Parallel Instances. The total number of solutions explored is set to 768k. The number of Parallel Instances is varied by with maximum n of 11, i.e. 1024 parallel instances. The best value after 2 -1 every 5 iterations is also recorded. We then compute the number of iterations required to reach a specific solution quality for different ACO architectures in Table 6 , expressed as proximity to the best-known optimal solution. For the particular problem and dataset, the best solution is the total cost of 2,701,367.58. There are six checkpoints of solution quality ranging from 99% to 99.9%. Although at first 1% gain might not seem significant, one must remember that global supply chain costs are measured in hundreds of millions, and even 1% savings do affect the bottom line. Empty fields (-) represent instances where the ACO was not able to converge to given solution quality. On all experiments, IAC was able to obtain solution quality only below 99.6%. In contrast, PA and PA with 5 ant local search were able to achieve above 99.9% solution quality with 512 and 1024 parallel instances. Furthermore, IAC did not see any significant benefit of adding more parallel instances for 99% and 99.25% checkpoints. (2,701,367.58) . Colour-coded from worse -in red, to the best -in green.

The number of Parallel Instances In contrast, PA does benefit from the increase in the number of parallel instances. For instance, PA can obtain the same solution quality in half the number of iterations at 99% checkpoint (scaling of 2x for sequential vs 1024 parallel instances). Scaling of 633.7x in case of 99.5% checkpoint for sequential counterpart. Similarly, PA with 5 ant sequential local search has the same dynamics, with scaling of 4x at 99% checkpoint compared to sequential and 140x at 99.6% checkpoint compared to 2 and 1024 parallel instances. One can also note that at increased solution quality and a little number of parallel instances, PA with 5 ant local search also offers improved efficiency in terms of total solutions explored. For example, at the 99.5% checkpoint with 2 parallel instances, PA takes 2590 iterations, while PA with 5 ant local search only requires 65 (decrease of 40x iterations or 8x total solutions explored). However, in most instances, PA without any local search is more efficient. 

To evaluate speed performance, we ran each given configuration and parallel architecture for 500 iterations or 10 minutes wall-clock time (whichever happens first) and recorded the total number of iterations and wall-clock time for three independent runs. Then, average wall-clock time per iteration was calculated. It is essential to measure the execution time correctly, just purely comparing computation per kernel/method may not show the real-life impact. For that reason, total time is measured from the start of the memory allocation to the freeing of the allocated memory, however it does not include the time required to load the dataset into memory. This allows us to estimate, with reasonable accuracy, what is the wall-clock time needed to run a specific architecture and configuration to converge to a given fitness quality. Although running each given architecture and configuration 10 times would produce more accurate convergence rate estimates, it would also require significantly more computation time. Furthermore, all vectorised implementations went through iterative profiling and optimisation process to obtain the fastest execution time. To the best of the authors' knowledge, all vectorised implementations have been fully optimised for the given hardware.

ACO implementation of IAC, PA and PAwV was implemented in C++ and multiple experiments of the configurations are shown in Table 7 . Intel C++ 18.0 with OpenMP 4.0 was used to compile the implementation. KMP 5 (an extension of OpenMP) config was varied based on total hardware core and logical core count (16c,2t = 32 OpenMP threads). Very similar results were obtained for both IAC double precision and PA double precision, with PA having around 5% overhead compared to IAC. In both instances, running 32 OpenMP threads offered around 24% speed reduction compared to 16 threads. Furthermore, PAwV with double precision vectorisation using AVX2 offered speed reduction of 26%, while scaling from 16 OpenMP threads to 32 offered almost no scaling at 256 parallel instances upwards. The nature of ACO pheromone sharing and probability calculations does not require double precision and therefore can be substituted with single-precision calculations. AVX2 offers 256-bit manipulations, therefore increasing theoretical throughput by a factor of 2, compared to double precision. 36% decrease in execution time was obtained, as not all parts of the code can take advantage of SIMD. Furthermore, doing 5 ant sequential local search within each parallel instance increases time linearly and produces little time savings in terms of solutions explored. The overall scaling factor at 1024 parallel instances compared to sequential execution at PAwV (single precision with AVX2 and 16c2t) is therefore 25.4x. 

Similar experiments were also conducted on the Xeon Phi hardware, Table 8 . Due to the poor convergence rate and search capability, the execution time for IAC was not measured. Xeon Phi differs from Hardware A with the ability to utilise up to 4 hyper-threads per core and AVX512 instruction set. Although Hardware B has 68 physical cores, for more straightforward comparison on base 2, only 64 were used in experiments. At 1024 parallel instances on double-precision PA, having 2 threads and 4 threads per core does offer speedup of 30% and 42% respectively, compared to 1 thread per core. Moving to the vectorised implementation of 256-bit AVX2, gains additional speedup of around 37% across all parallel instances, however, did not benefit from 4 hyper-threads. Furthermore, exploiting the AVX512 instruction set offers a further 24% speedup compared to AVX2. In this configuration having 4 hyper threads per core worsens the speed performance (3.644 seconds vs 3 seconds). Like Hardware A, PAwV was explored with single precision and offered near-perfect scaling on 1024 parallel instances with 4 hyper-threads per core, or 40% overall speed improvement compared to PAwV with double precision (3 seconds vs 1.804 seconds). Alike Hardware A, having 5 sequential local ants does not provide any time savings and time increases linearly. The overall scaling factor at 1024 parallel instances compared to sequential execution at PAwV (single precision with AVX512 and 64c4t) is therefore 148x. 

A further set of experiments were also conducted for GPU, Table 9 . The implementation with no vectorisation (Blocks x1), uses 1 thread per CUDA block to compute one solution, therefore 1024 parallel instances require 1024 blocks. Similarly, for (Blocks x32), 32 threads are used per block, each thread computing its own solution independently. For parallel instances of 32, only 1 block would be used with 32 threads. The implementation of no vectorisation utilises no shared memory; however, all static problem metadata is stored as textures. A single kernel is launched, and the best solution across all parallel instances is returned. Vectorized version implements architecture described in [68] , storing the route choice matrix in shared memory and utilising local warp reduction for sum and max operations. Each thread-block builds its solution, while the extra 32 threads assist with the reduction operations, memory copies and fitness evaluation. Table 9 shows a comparison between the two implementations. Implementation without vectorisation performs on average two times slower compared to the vectorised version. Furthermore, 64 threads per block (Blocks x64) performs slower than 32 threads per block (Block x32). Next, scaling across multiple GPUs were explored. Each device takes a proportion of 1024 instances with unique seed values and after each iteration, the best overall solution is reduced. In the case of 2 GPUs and 1024 parallel instances, each device will compute 512 parallel instances concurrently. Scaling across 2 (2x) and 4 GPUs (4x) did not provide any significant speedup (only 10%). This is due to the fact that each iteration consumes at least 50 seconds and scaling across multiple GPUs adds almost no overhead. The maximum number of parallel instances might need to be increased to fully utilise all 4 GPUs to the point where all Streaming Multiprocessors (SMs) are saturated and increasing block count increases the computation time linearly. GPU implementation is, therefore, one magnitude of order slower than that of CPU. However, this could be explained by the nature of the problem and not be specific to ACO architecture, as there have been a lot of success on GPUs solving simple, low memory footprint TSP instances [57] [68] [70] . However, the problem addressed in this paper requires a lot of random global memory access to check for all restrictions such as order limits, capacity constraints and weight limits, which are too big to be stored on the shared memory. 

If both convergence rate of the architecture and the speed of the hardware is considered, an estimate can be made on what would be the average wall-clock time to converge to specific solution quality. The fastest configuration for both Hardware A (Table 7) and Hardware B (Table 8 ) was chosen and then multiplied by the number of iterations required to reach a specific solution quality (   Table 6 ) to obtain an estimate of the compute time required (Table 10) In addition, best computation time to achieve specific solution quality was also compared in Figure 7 , where the estimated best computation time required (in logarithmic) is plotted against three tested architectures across various solution quality checkpoints. Figure 7 clearly shows that GPU results (Hardware C) were considerably slower and therefore, authors conclude that GPUs are not suitable for the supply chain problem solved. 

Solution quality Estimated time required (in seconds, log) 

In addition to the real-world supply chain problem, a single TSP instance with 318 cities (lin318) is selected for comparison. The lin318 instance is small enough such that all experiments can be computed quickly but large enough to see measurable differences between hardware architectures explored. Like in the supply chain problem, solution quality checkpoints against optimal fitness value of 42029 were recorded during the convergence process. Moreover, just like in supply chain problem, PA outperformed IAC architecture for solving lin318. The lin318 computation time was plotted against various hardware solutions and solution quality checkpoints in Figure 8 . When solving the lin318 TSP instance, Hardware A performs faster than Hardware B for solution quality between 99.0 and 99.6% and slower for higher solution quality, similar to the supply chain problem results in Figure 7 . Although Hardware C -GPU performed magnitudes slower in supply chain problem, for the TSP instance it was able to converge faster than Hardware A and Hardware B. Therefore, authors can confirm the findings of [57] [68] [70] , that suggest that GPUs offer speedup over CPU counterpart when routing simple TSPs. However, authors also acknowledge that these dynamics do not apply for a more complex real-world routing problem where GPU is magnitudes slower than CPU counterparts (Hardware A or Hardware B) due to the additional meta-data required to be stored during solution creation.

Nature-inspired meta-heuristic algorithms such as Ant Colony Optimization (ACO) have been successfully applied to multiple different optimisation problems. Most work focuses on the Travelling Salesman Problem (TSP). While TSPs are a good benchmark for new idea comparison, the dynamics of the proposed algorithms for benchmarks do not always match to a real-world performance where the problem has more constraints (more metadata during solution creation). Furthermore, speed and fitness performance comparisons are not always completely fair when compared to a sequential implementation. This work explored the dynamics of different ACO architectures applied to benchmark and real-world problem. The experimental results demonstrate that the results obtained from the TSP benchmarks cannot be generalised to the real-world applications, especially in terms of hardware performance and usage. Therefore, our findings demonstrate that in order to achieve the generalisable conclusions, the experimental work has to be completed on both: standard benchmarks and real-world applications. Furthermore, our work solves a real-world outbound supply chain network optimisation problem and compares two different ACO architectures -Independent Ant Colonies (IAC) and Parallel Ants (PA). It was concluded that PA outperformed IAC in all instances, as IAC failed to find any better solution than 99.5% of optimal. In comparison, PA was able to find a near-optimal solution (99.9%) in fewer iterations due to effective pheromone sharing across ants after each iteration. Furthermore, PA shows that it consistently finds a better solution with the same number of iterations as the number of parallel instances increase. Moreover, a detailed speed performance was measured for three different hardware architectures -16 core 32 thread workstation CPU, 68 core server-grade Xeon Phi and general-purpose Nvidia GPUs. Results showed that although GPUs can scale when solving simple TSP (as confirmed by multiple other studies), those scaling dynamics do not transfer to more complex the real-world problem. Memory access footprint required to check capacity limits and weight constraints did not fit on small shared memory on GPU. Thus, it performed 29 times slower than the other two hardware solutions even when running 4 GPUs in parallel. Therefore, authors consider this finding a new knowledge with surprise value. When compared to a real-life impact on the time required to reach a specific solution quality, both CPU and Xeon Phi optimised-vectorised implementations showed comparable speed performance; with CPU taking the lead with lower Parallel Instances count due to the much higher clock frequency. At near-optimal solution (99.75%+) and 1024 parallel instances, Xeon Phi was able to take full advantage of AVX512 instruction set and outperformed CPU in terms of speed. Therefore, compared to an equivalent sequential implementation at 1024 parallel instances, CPU was able to scale 25.4x while Xeon Phi achieved a speedup of 148x. Since PA fitness performance increases as the number of parallel instances increase, it would be worth looking into scaling above 1024 instances using either cluster of CPUs or clusters of Xeon Phi's, which will be part of the future work. Furthermore, Field Programmable Gate Arrays (FPGAs) might have the potential to take advantage of highly vectorised ACO, which is another area of possible future research.

646% 98.701% 98.740% 98.713% 98.813% 98.825% 98.857% 98.859% 98.881% 98.931% 98.923% 20 98.921% 98.935% 98.973% 98.987% 98.980% 99.063% 99.053% 99.082% 99.102% 99.133% 99.150% 40 99.165% 99.265% 99.315% 99.300% 99.343% 99.355% 99.366% 99.413% 99.410% 99.427% 99.443% 60 99.354% 99.413% 99.466% 99.503% 99.530% 99.536% 99.541% 99.562% 99.573% 99.592% 99.598% 80 99.438% 99.459% 99.547% 99.547% 99.585% 99.585% 99.582% 99.630% 99.638% 99.660% 99.667% 100 99.444% 99.459% 99.548% 99.559% 99.589% 99.592% 99.584% 99.646% 99.641% 99.672% 99.674% 200 99.452% 99.461% 99.551% 99.569% 99.601% 99.605% 99.599% 99.724% 99.717% 99.846% 99.844% 300 99.452% 99.461% 99.558% 99.574% 99.615% 99.615% 99.606% 99

Tactical supply chain planning models with inherent flexibility: definition and review

An ant colony system for responsive dynamic vehicle routing

Multi-satellite control resource scheduling based on ant colony optimization

A two-phase genetic algorithm for incorporating environmental considerations with production, inventory and routing decisions in supply chain networks

Using multi-objective genetic algorithm for partner selection in green supply chain problems

A green home health care supply chain: New modified simulated annealing algorithms

A Meta-Heuristic Algorithm Based on Simulated Annealing for Designing Multi-Objective Supply Chain Systems

Industrial & Systems Engineering Conference (ISEC)

An ant colony system empowered variable neighborhood search algorithm for the vehicle routing problem with simultaneous pickup and delivery

Ant Colony Algorithm for Routing Alternate Fuel Vehicles in Multi-depot Vehicle Routing Problem

Resilient food supply chain design: Modelling framework and metaheuristic solution approach

Development of an ant colony optimisation-based heuristic for a location-routing problem in a two-stage supply chain

A Comparison of ACO, GA and SA for Solving the TSP Problem

HOPNET: A hybrid ant colony optimization routing algorithm for mobile ad hoc network

A novel hybrid approach based on Particle Swarm Optimization and Ant Colony Algorithm to forecast energy demand of Turkey

A revised ant algorithm for solving location-allocation problem with risky demand in a multi-echelon supply chain network

Ant Colony Optimization For Split Delivery Inventory Routing Problem

Designing closed-loop supply chains with nonlinear dimensioning factors using ant colony optimization

MMAS on GPU for Large TSP Instances

A GPU-based parallel MAX-MIN Ant System algorithm with grouped roulette wheel selection

A Parallel Implementation of Ant Colony Optimization

Metaheuristic algorithms and probabilistic behaviour: a comprehensive analysis of Ant Colony Optimization and its variants

A parallel cooperative hybrid method based on ant colony optimization and 3-Opt algorithm for solving traveling salesman problem

Parallel Performance of an Ant Colony Optimization Algorithm for TSP

A Survey on GPU-Based Implementation of Swarm Intelligence Algorithms

First results of performance comparisons on many-core processors in solving QAP with ACO

Parallel ant colony optimization for resource constrained job scheduling

Comparative evaluation of platforms for parallel Ant Colony Optimization

RMACO :a randomly matched parallel ant colony optimization

Parallel ant colony optimization on multi-core SIMD CPUs

Parallel Ant Colony Optimization for the HP Protein Folding Problem

A CUDA based Solution to the Multidimensional Knapsack Problem Using the Ant Colony Optimization

CPU Versus GPU Parallelization of an Ant Colony Optimization for the Longest Common Subsequence Problem

High-throughput Ant Colony Optimization on graphics processing units

Implementing a GPU-based parallel MAX-MIN Ant System

Accelerating ant colony optimization-based edge detection on the GPU using CUDA

Parallel Ant Colony Optimization for Flow Shop Scheduling Subject to Limited Machine Availability

Using a coprocessor to solve the Ant Colony Optimization algorithm

Efficient exploitation of the Xeon Phi architecture for the Ant Colony Optimization (ACO) metaheuristic

A highly parallelized and vectorized implementation of Max-Min Ant System on Intel® Xeon Phi TM

An Effective Parallelism Topology in Ant Colony Optimization algorithm for Medical Image Edge Detection with Critical Path Methodology (PACO-CPM)

Multi-threading based implementation of Ant-Colony Optimization algorithm for image edge detection

Applying ACO to Large Scale TSP Instances

Effective heuristics for ant colony optimization to handle large-scale problems

Performance analysis and tuning for parallelization of ant colony optimization by using openmp

Parallel ant colony optimization for the determination of a point heat source position in a 2-D domain

A Parallel Ant Colony Optimization for the Maximum-Weight Clique Problem

Evaluation of Parallelism in Ant Colony Optimization Method for Numerical Solution of Optimal Control Problems

Generic Techniques in General Purpose Gpu Programming With Applications To Ant Colony and Image Processing Algorithms

Accelerating Ant Colony Optimization for the Vertex Coloring Problem on the GPU

Optimizing PolyACO Training with GPU-Based Parallelization

GPU implementation of ant colony optimization-based band selections for hyperspectral data classification

A GPU-based Parallel Ant Colony Algorithm for Scientific Workflow Scheduling

Transit stop inspection and maintenance scheduling: A GPU accelerated metaheuristics approach

Research on Solving Travelling Salesman Problem using Rank Based Ant System on GPU

Multi Colony Ant System based Solution to Travelling Salesman Problem using OpenCL

Query Optimization using Modified Ant Colony Algorithm

Accelerating ant colony optimisation for the travelling salesman problem on the GPU

Dynamic strategy based parallel ant colony optimization on GPUs for TSPs

Machine Learning, Optimization, and Big Data

ACO-PSO optimization for solving TSP problem with GPU acceleration

Dynamic load balancing on heterogeneous clusters for parallel ant colony optimization

Enhancing data parallelism for Ant Colony Optimization on GPUs

The GPU-based parallel Ant Colony System

Parallel Programming in OpenMP

Intel Xeon Phi Processor High Performance Programming: Knights Landing Edition

Supply Chain Logistics Problem Dataset

Ant colony system: a cooperative learning approach to the traveling salesman problem

High-throughput Ant Colony Optimization on graphics processing units

Composite goal methods for transportation network optimization

GACO: A GPU-based High Performance Parallel Multi-ant Colony Optimization Algorithm

Authors would like to thank Intel Corporation for donating the Xeon Phi hardware.

Authors would like to thank Intel Corporation for donating the Xeon Phi hardware.

 Standard TSP instances are not suitable for generalized conclusions  Although good for TSPs, GPUs are not well suited for explored supply chain problem  25.4x speed-up was achieved on CPU compared to its sequential counterpart  148x speed-up was achieved on Xeon Phi compared to its sequential counterpart