Parallel Computing Systems for the Optimal Selection of Process Equipment

arXiv:2611.16290v1 [cs.DC] 12 Nov 2026

A. B. BORISENKO

PARALLEL COMPUTING SYSTEMS FOR THE OPTIMAL SELECTION OF PROCESS EQUIPMENT

Monograph English edition of the author’s manuscript

Tambov 2026

UDC 004.272.2:519.854:66.011

The Russian edition of this monograph was recommended for publication by the Scientific and Technical Council of Tambov State Technical University

Reviewers of the Russian edition: E. N. Tugolukov, Doctor of Engineering Sciences, Professor, Department of Engineering and Technology of Nanoproduct Manufacturing, Tambov State Technical University; A. D. Obukhov, Doctor of Engineering Sciences, Associate Professor, Department of Automated Decision Support Systems, Tambov State Technical University; D. S. Solovjev, Doctor of Engineering Sciences, Associate Professor, Department of Mathematical Modeling and Information Technologies, Derzhavin Tambov State University

Borisenko, A. B. Parallel Computing Systems for the Optimal Selection of Process Equipment : monograph / A. B. Borisenko. — Tambov, 2026. — 99 p. The monograph is devoted to the development and study of parallel algorithms for the optimal selection of equipment for chemical-engineering systems of multiproduct batch plants, which is a large-scale NP-hard discrete optimization problem. It presents a mathematical model of the operation of a chemicalengineering system; sequential and parallel variants of the branch-and-bound method for shared-memory and distributed-memory systems (OpenMP, MPI) and for graphics processing units (CUDA); a hybrid metaheuristic algorithm combining Ant Colony Optimization and Simulated Annealing; and methods for tuning its parameters on one or several GPUs. Results of computational experiments on multi-core processors, computing clusters, and GPUs of various architectures are reported. The book is intended for researchers, engineers, and software developers working in parallel and high-performance computing, as well as for doctoral and senior undergraduate students in related fields. This English edition is the author’s translation of the Russian-language monograph Primenenie parallel’nykh vychislitel’nykh sistem dlya optimal’nogo vybora tekhnologicheskogo oborudovaniya (Tambov, 2026), with minor corrections.

© A. B. Borisenko, 2026

Contents Abbreviations

5

Introduction

6

1 Optimal Equipment Selection for Chemical-Engineering Systems: Problem Statement 8 1.1 Mathematical Formulation of the Problem . . . . . . . . . . . . . . . . . . . . . . 9 1.2 Mathematical Model of CES Operation . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.1 Input Data and Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.2 Batch Size and the Compatibility Constraint . . . . . . . . . . . . . . . . 12 1.2.3 Operation Durations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.2.4 Schedule of a System Cycle . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2.5 Production Time and the Processing Time Constraint . . . . . . . . . . . 17 1.2.6 Optimality Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.2.7 Computational Aspects of the Model . . . . . . . . . . . . . . . . . . . . . 18 1.3 Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Sequential Branch-and-Bound Algorithm 2.1 Basic Version of the Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Execution Times of Individual Parts of the Algorithm . . . . . . . . . . . . . . . 2.3 Optimization of the Sequential Algorithm . . . . . . . . . . . . . . . . . . . . . .

21 21 23 24

3 Parallel Branch-and-Bound for Multi-Core and Cluster Systems 3.1 Methods for Parallelizing Search Tree Traversal . . . . . . . . . . . . . . . . . . . 3.1.1 Features of Parallel Search Tree Traversal . . . . . . . . . . . . . . . . . . 3.1.2 Static and Dynamic Tree Partitioning . . . . . . . . . . . . . . . . . . . . 3.1.3 Classification of Dynamic Load Balancing Methods . . . . . . . . . . . . . 3.1.4 Algorithmic Aspects of Dynamic Load Balancing . . . . . . . . . . . . . . 3.1.5 Speedup Anomalies and the Search Overhead Factor . . . . . . . . . . . . 3.1.6 The Assigned Subtrees and Allocated Subtrees Methods . . . . . . . . . . 3.1.7 Parallelizing Branch-and-Bound on GPUs . . . . . . . . . . . . . . . . . . 3.1.8 Methods Used in This Book . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Shared-Memory Systems (OpenMP) . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Program Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Results of Computational Experiments . . . . . . . . . . . . . . . . . . . . 3.3 Distributed-Memory Systems (MPI) . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Master Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Worker Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Results of Computational Experiments . . . . . . . . . . . . . . . . . . . .

25 25 25 26 26 27 28 29 29 30 30 30 31 32 32 33 34

3

4 Branch-and-Bound on Graphics Processing Units 4.1 General Principles and Features of GPU Programming . . . . . . . . . . . . . . . 4.2 Development of the Parallel Program . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Optimization of the Program Code . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Results of Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . .

37 38 39 44 47

5 Parallel Hybrid Metaheuristic Algorithm on GPUs 5.1 A Brief Overview of Metaheuristic Algorithms for Combinatorial Optimization . 5.2 Simulated Annealing (SA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Ant Colony Optimization (ACO) . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 GPU Implementation of the Hybrid Algorithm (ACO+SA) . . . . . . . . . . . . 5.5 Results of Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . .

51 51 53 55 58 63

6 Tuning the Parameters of the Ant Colony Algorithm 6.1 Approaches to Tuning the Parameters of Metaheuristic Algorithms . . . . . . . . 6.2 Parameters of the Ant Colony Algorithm . . . . . . . . . . . . . . . . . . . . . . 6.3 Proposed Tuning Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Results of Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Experiments on GPUs of Different Architectures and on Multiple GPUs . . . . . 6.5.1 Organization of Computations on Multiple GPUs . . . . . . . . . . . . . . 6.5.2 Experimental Platforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.3 Comparison of GPU Architectures . . . . . . . . . . . . . . . . . . . . . . 6.5.4 Scaling the Parallel Program to Multiple GPUs . . . . . . . . . . . . . . .

68 68 70 71 73 77 78 79 80 80

Conclusion

84

Bibliography

86

4

Abbreviations ACO API B&B BFS CC CES CO CPU CSP CUDA DFS DOE DRAM DSP FP32 FPGA GPGPU GPU HLSL LIFO MAC MILP MINLP MIR MPI NP NVCC OFAT OpenCL OpenMP RAM SA SIMD SIMT SM SMP STL

Ant Colony Optimization Application Programming Interface Branch-and-Bound Breadth-First Search Compute Capability Chemical-Engineering System Combinatorial Optimization Central Processing Unit Constraint Satisfaction Problem Compute Unified Device Architecture Depth-First Search Design of Experiments Dynamic Random-Access Memory Digital Signal Processor 32-bit (single-precision) Floating Point Field-Programmable Gate Array General-Purpose computing on Graphics Processing Units Graphics Processing Unit High-Level Shading Language Last-In, First-Out Multiple Ant Colonies Mixed-Integer Linear Programming Mixed-Integer Nonlinear Programming Multiple Independent Runs Message Passing Interface Nondeterministic Polynomial (time) NVIDIA CUDA Compiler One-Factor-At-a-Time Open Computing Language Open Multi-Processing Random-Access Memory Simulated Annealing Single Instruction, Multiple Data Single Instruction, Multiple Threads Streaming Multiprocessor Symmetric Multiprocessing Standard Template Library

5

Introduction NP-hard combinatorial optimization (CO) problems are among the most difficult computational problems to solve in practice. A CO problem is an optimization problem in which a finite set of objects is given together with an objective function that assigns a cost to each object. The goal is to find an object with the minimum (or maximum) value of the objective function. The objects are typically integers, subsets of a set of elements, permutations of a set of elements, and so on. Well-known examples of CO problems include the Travelling Salesman Problem [1], the Assignment Problem [2], the Knapsack Problem [3], and various scheduling problems [4]. A characteristic feature of NP-hard problems is that, when exact (deterministic) algorithms such as exhaustive enumeration are used, the time required to find the optimal solution grows exponentially with the size of the search space; this phenomenon is known as combinatorial explosion [5]. The main ways of increasing computational speed are the following: using more advanced algorithms; tuning the parameters of the algorithms employed; compiler optimization of the programs; optimization with the help of profilers; and using more powerful hardware, including the development of parallel computing technologies. Of these, the strongest influence on computational speed comes from the choice of the solution algorithm together with an optimal setting of its parameters, and from the hardware used. All known algorithms for NP-hard problems, including CO problems, can be roughly divided into three categories: exact [6], heuristic [7], and metaheuristic [8] algorithms. Exact algorithms include, for example, exhaustive search [9], the branch-and-bound method [10], the Hungarian algorithm [11], and dynamic programming [12]. Heuristic approaches rely on plausible but not mathematically justified assumptions about the properties of the optimal solution [13]. Such algorithms usually yield practically good solutions, but they do not guarantee that the best solution will be found. Moreover, heuristics are built on “common sense” and therefore depend strongly on the particular problem. Metaheuristics combine basic heuristic methods within higher-level algorithmic frameworks aimed at an efficient exploration of the search space [14–16]. Among metaheuristic algorithms are Tabu Search [17], Simulated Annealing [18], Ant Colony Optimization [19], Genetic Algorithms [20], Particle Swarm Optimization [21], Artificial Bee Colony Optimization [22], and others. Their advantage over classical heuristics based on local search is that they can explore a larger part of the search space in order to find a near-optimal solution, whereas local search methods stop as soon as a local optimum in the chosen neighborhood has been reached. It should also be noted that hybridization of metaheuristics is currently one of the most promising ways of designing new, highly efficient algorithms [23, 24]. A distinctive feature of metaheuristic algorithms is that they contain a large number of free parameters. The efficiency of an algorithm depends strongly on the values of these parameters; however, as a rule, there are no universal formal guidelines for choosing them on the basis of the properties of the optimization problem at hand [25]. The development of methods for tuning metaheuristics is therefore a topical research problem [26, 27]. Over the past decades, parallelization of algorithms and the use of high-performance computing systems have become one of the most promising ways of improving algorithm performance [28–31]. Such systems include shared-memory systems, distributed-memory systems [32], and heteroge6

neous systems. Shared-memory systems include computers with one or more multi-core processors. Distributed-memory systems consist of several compute nodes, each with its own main memory, connected into a computing cluster by communication links (e.g., Ethernet or InfiniBand). Most modern high-performance computing systems have a heterogeneous architecture comprising a general-purpose multi-core Central Processing Unit (CPU) and a set of various accelerators [33– 35]. The performance of heterogeneous systems is increased not only by increasing the number of processors and cores, but also by incorporating (either integrated or discrete) specialized accelerators for parallel computing, such as Graphics Processing Units (GPUs), Digital Signal Processors (DSPs), and Field-Programmable Gate Arrays (FPGAs) [36]. Traditional approaches to developing parallel programs for CPU-based systems are not sufficient here. Efficient implementation of algorithms on heterogeneous systems is possible only if the specific features and limitations of their architecture are taken into account, including the interaction between the CPU and the accelerator, the memory hierarchy and the rules for using it efficiently, and the orientation towards massively parallel data processing. Thus, because of the high computational complexity of practically important optimization problems and the rapidly growing availability of high-performance computing systems of various kinds, the development of parallel metaheuristic algorithms is becoming increasingly important. This book is devoted to the development and study of parallel discrete optimization algorithms, the tuning of their parameters, and their efficient implementation on various parallel computing systems, including distributed-memory systems, shared-memory systems, and heterogeneous systems equipped with GPUs. As a test problem, we use the optimal selection of equipment for chemical-engineering systems (CES) [37]: among all available alternatives, one has to find the combination of equipment units for which the objective function (the total capital cost of the equipment) attains its minimum. The book is organized as follows. Chapter 1 states the problem. Chapter 2 presents a sequential algorithm for solving it based on the branch-and-bound (B&B) method. Chapter 3 presents parallel implementations of the algorithm for shared-memory and distributed-memory systems. Chapter 4 is devoted to the implementation of branch-and-bound on heterogeneous systems equipped with GPUs. Chapter 5 develops a parallel hybrid metaheuristic algorithm for the problem and its implementation on heterogeneous systems. Chapter 6 deals with the selection and optimization of the parameters of the hybrid metaheuristic algorithm.

7

Chapter 1

Optimal Equipment Selection for Chemical-Engineering Systems: Problem Statement Selecting the equipment of a chemical-engineering system (CES) is one of the main problems that arise both in designing new [38] and in retrofitting existing [39] multiproduct batch plants. Such plants are characterized by relatively small production volumes, a wide range of products, high product prices, equipment operating in batch or semi-continuous mode and, accordingly, production in batches, the ability to switch quickly from one product to another using the same set of equipment (high flexibility), and a large number of processing stages. Typical examples of multiproduct batch plants are the production of chemical reagents, pharmaceuticals, cosmetics, dyes, and intermediates [40, 41]. At each processing stage, the materials are processed for a certain time and then transferred to the next stage. Batches of each product must pass successively through the equipment of all production stages, and the units of each stage must be able to receive and process a batch completely. For example, a batch of material from a unit with a volume of 10 m3 cannot be transferred to a unit with a working volume of 0.1 m3 . If such a situation arises, the CES is not operable. The operating cycle of a unit is an ordered sequence of main and auxiliary operations of finite duration. The most common of them are loading, the main processing operations, unloading, and cleaning with preparation for the next cycle. In addition, a unit may be idle during operation. The processing time at a stage depends on the main size (working volume or working surface area) of the unit in which the processing takes place, because a change in the unit size changes the amount of material processed in it. Therefore, whenever the main sizes of the equipment change, the production time T has to be recalculated. The plant being designed must produce the required amount of products Q within a given time horizon Tmax . A solution of the multiproduct plant design problem is a CES variant, i.e., the number of main units at each stage of the CES and the main size (working volume or working surface area) of each unit, together with the operating parameters (product batch sizes and production periods) that ensure the required production rate for all products. The main size is usually chosen from a set of standard values. One has to find the optimal combination of the number and main sizes of units at the stages of the CES among the set of possible variants, together with the operating parameters of a system with the given production rate; the optimality criterion may be, for example, the total capital cost of the equipment [42–44]. The solution of this problem largely determines the quality of design decisions on equipment layout and production scheduling. In this book, we consider the selection of the main equipment of the CES stages. The selection of auxiliary equipment (measuring tanks, receivers, pumps, etc.) is a separate problem, which is solved using the results 8

of the main equipment selection [45, 46]. The problem under consideration belongs to the class of Mixed-Integer Nonlinear Programming (MINLP) problems [47–49]. In a number of works, it is reduced to a Mixed-Integer Linear Programming (MILP) problem [50, 51]. For each operable CES variant, the operating parameters of the system are determined uniquely. Therefore, the search for the optimal CES variant can be regarded as a discrete optimization problem. Finding the global optimum by exhaustive enumeration of all variants is impractical because of the large size of the problem. As the size of the search space grows (for this problem, as the number of available standard equipment sizes increases), the number of combinations to be examined (all possible system variants E) grows exponentially because of the combinatorial explosion, so that even the most powerful computers may require an unacceptable amount of time to find a solution. Various approaches to this problem have been proposed in the literature, for example, Monte Carlo methods, genetic algorithms, heuristic methods, and the reduction of the original problem to a nonlinear or linear programming problem [38, 47, 48, 52, 53]. Heuristic approaches aim at finding acceptable (suboptimal) solutions, whereas mathematical programming methods guarantee exact results. However, as the problem size grows, exact methods face the combinatorial explosion, which makes them inefficient. In this book, we propose parallel algorithms for the equipment selection problem for CES of multiproduct batch plants and study their applicability to finding the optimal solution. The chapter is organized as follows. Section 1.1 gives the mathematical formulation of the optimal CES equipment selection problem and describes the representation of the search space as a tree of variants. Section 1.2 considers the mathematical model of CES operation, which determines the form of the constraints on the operability and on the production time of the system as well as of the optimality criterion, and analyzes the computational properties of the model that are essential for designing sequential and parallel algorithms. Section 1.3 describes the test case used in the computational experiments throughout the subsequent chapters.

1.1

Mathematical Formulation of the Problem

A CES of a multiproduct batch plant is a set of processing stages equipped with reactors, tanks, filters, dryers, and other process equipment, and intended to carry out the production steps for a given range of products. Assuming that one main unit can be installed at each processing stage of the CES, the problem can be formulated as follows [54, 55]. A CES is a sequence of I processing stages. Each stage i of the system can be equipped with a main unit chosen from a finite set Xi containing Ji elements. The set Xi comprises all possible main sizes of the units for stage i: Xi = {xi,j }, i = 1, I, j = 1, Ji , where xi,j is the main size (working volume or working surface area) of the jth unit that can be installed at stage i of the CES. Q Each CES variant Ωe , e = 1, E, where E = Ii=1 Ji , is an ordered set of main sizes of the units chosen from the corresponding sets. Figure 1.1 shows an example of a simple CES consisting of 4 stages (I = 4). For each processing stage of the CES, two standard unit sizes are available, i.e., J1 = J2 = J3 = J4 = 2. Here X1 = {x1,1 , x1,2 }, X2 = {x2,1 , x2,2 }, X3 = {x3,1 , x3,2 }, X4 = {x4,1 , x4,2 }. The total number of possible CES variants is E = J1 · J2 · J3 · J4 = 24 = 16. The order in which product batches are processed by the main units of the CES stages is fixed (it is determined in advance when the structure of the material flows is designed); for example, the variants {x1,1 , x2,1 , x3,2 , x4,2 } and {x1,2 , x2,1 , x3,1 , x4,1 } are admissible, whereas the variants {x3,1 , x2,1 , x4,1 , x1,2 } and {x2,2 , x4,2 , x3,1 , x1,1 } are obviously inadmissible. Every operable CES variant must satisfy the conditions for the joint operation of the equipment of its stages, S(Ωe ) = 0 (the compatibility constraint). A CES variant is operable if the units of 9

Figure 1.1: Example of a simple chemical-engineering system every one of its processing stages, without exception, are able to process the batch volume of the product, taking into account the maximum and minimum filling degrees of the units. In addition, an operable CES variant must ensure the production of the given product range and amounts within a given period of time, i.e., T (Ωe ) ≤ Tmax (the processing time constraint), where Tmax is the given time horizon. Thus, the CES equipment selection problem can be formulated as follows: on the set of available process equipment (1.2)–(1.3), find a CES variant Ω∗ ∈ Ωe , e = 1, E, for which the optimality criterion, the capital cost Cost(Ωe ) (1.1), attains its minimum and the compatibility constraint (1.4) and the processing time constraint (1.5) are satisfied: Ω∗ = argmin Cost(Ωe ), e = 1, E

(1.1)

Ωe = {x1,j1 , x2,j2 , . . . , xI,jI |ji = 1, Ji , i = 1, I}, e = 1, E

(1.2)

xi,j ∈ Xi , i = 1, I, j = 1, Ji

(1.3)

S(Ωe ) = 0, e = 1, E

(1.4)

T (Ωe ) ≤ Tmax , e = 1, E

(1.5)

The search space for the CES shown in Fig. 1.1 is a tree consisting of I levels (see Fig. 1.2).

Figure 1.2: Search space 10

Each level i of the tree corresponds to a stage of the CES. Each edge xi,j of the tree corresponds to the selected unit variant j for stage i (one selected element of the set Xi = {xi,j }, i = 1, I, j = Q 1, Ji ). Each internal node ni,k , i = 1, I, k = 1, Ki , Ki = il=1 (Jl ) of the tree corresponds to a fragment of the CES consisting of stages 1 to i. Note that the problem under consideration is NP-hard [48, 56, 57]. The form of the functions S(Ωe ) and T (Ωe ), as well as of the optimality criterion, is determined by the mathematical model of CES operation discussed in Section 1.2.

1.2

Mathematical Model of CES Operation

In the problem formulation (1.1)–(1.5), the compatibility and processing time constraints are written in terms of the functions S(Ωe ) and T (Ωe ), whose specific form is determined by the mathematical model of system operation [58]. The model is presented here to the extent needed to understand the computational side of the problem: which quantities are computed when a single CES variant is evaluated, which of these computations can be performed for incomplete variants (internal nodes of the search tree), and which data remain unchanged throughout the search. The process engineering justification of the model relations and the procedure for preparing the input data from process specifications are beyond the scope of this book; they are described in detail in [58]. Models of the operation of multiproduct batch plants are widely represented in the literature [59, 60]. In the most common formulation, each stage is characterized by a size factor (a counterpart of the material index used below) and a batch processing time, the production rate of the system is determined by the limiting cycle time, and the time constraint is written as the sum, over all products, of the number of batches multiplied by the corresponding limiting cycle time [61, 62]. Within this formulation, the literature considers parallel units operating in phase or out of phase [63], intermediate storage tanks and discrete sets of standard equipment sizes [62], as well as units duplicated in series [64]. The model used in this book [58] is detailed down to the individual operations of the unit cycles: it explicitly accounts for the overlap of loading and unloading with the processing operations, for the splitting and merging of batches at particular stages, and for the resulting “filled idle times” of the equipment. The most complete description of the model and of the methods for solving its relations is given in [65–67]. The notation is consistent with Section 1.1: i = 1, I is the index of a processing stage of the CES, j = 1, Ji is the index of a standard unit size, and Ωe is a CES variant. The main size of the unit selected for stage i in the variant Ωe = {x1,j1 , . . . , xI,jI } will be denoted by x̂i = xi,ji .

1.2.1

Input Data and Assumptions

The CES manufactures a range of P products, p = 1, P . For each product, the production amount Qp (t) is given; the total time available for producing the whole product range is Tmax . Product p is processed at a subset of processing stages Ip ⊆ {1, . . . , I}. The processing order is defined by the process route, a directed graph Γp whose vertices are the stages in Ip and whose arcs correspond to the transfer of materials from stage to stage. Let Pp,i denote the set of stages from which materials are fed to stage i when product p is produced; stages with Pp,i = ∅ will be called initial stages. Dye production typically has routes in the form of a converging tree: several independent branches in which raw materials and intermediates are prepared merge at the stages of chemical transformation, and the route ends with a single final stage. For example, in the test case (Section 1.3), when the first product is manufactured, stage 110 (azo coupling) receives materials from stages 103 and 106, and stage 145 (drying) is the final stage for all three products. According to the type of the main unit, the set of stages is partitioned into four disjoint classes: A, batch vessels (reactors and tanks equipped with agitators); B, filters and dryers, in which the processing time is determined by the amount of material processed and by the working

11

surface area (or volume) of the unit; C, chamber and plate-and-frame filter presses whose target product is the solid phase of the suspension; D, batch vacuum dryers. Stage i may be equipped with νi ≥ 1 identical units of the same standard size x̂i ; the number νi is part of the input data and is not an optimization variable. For νi > 1, the units operate either out of phase (each unit processes a whole batch in turn) or in phase (a batch is divided equally among the units). The model is based on the following assumptions. 1. At any moment, the CES manufactures only one product; the products of the range are manufactured one after another. 2. The operating cycle of a unit at any stage is an ordered sequence of operations of six types: loading, “filled idle time” during loading, physicochemical transformations (hereafter referred to as processing), unloading, “filled idle time” during unloading, and cleaning (preparation for the next cycle). A “filled idle time” is a period during which the unit stays partially loaded or partially unloaded, waiting for the next portion of material or for the receiving unit to become ready. 3. Loading and unloading are either performed as separate operations or overlap in time with processing (for example, when a suspension is fed continuously to a filter). 4. The duration of each operation is either constant or a known function of the amount of material processed and of the main size of the unit; it does not depend on the batch number or on the cycle number. 5. At particular stages, a batch may be divided into equal portions processed in succession, or several batches may be merged and processed together. Note that the FPAKM automatic chamber filter presses installed at the filtration stages of the test CES (stages 108, 111, and 134) belong to class B despite their name: they operate cyclically like ordinary filters (suspension feeding, filtration, cake discharge, filter cloth regeneration), and the filtration time is determined by their working surface area and the batch size. In the model, loading and unloading of these units overlap with filtration, while the time for cake discharge and cloth regeneration is accounted for by the fraction 1 − hp,i of auxiliary operations (see (1.14)). The test case contains no stages of classes C and D. The way in which a batch of product p is processed at stage i is specified by the batch size ratio rp,i : for rp,i = 1, the batch is processed as a whole; for rp,i = γ > 1, it is divided into γ equal portions; for rp,i = 1/γ, γ > 1, γ batches are merged at the stage. The fraction of a batch that resides in one unit of the stage at a time is given by the coefficient (

up,i =

1/rp,i , νi = 1 or out-of-phase operation, 1/(νi rp,i ), in-phase operation.

(1.6)

The process specification for product p provides the following characteristics for each stage i ∈ Ip : vp,i and µp,i , the volumetric (m3 /t) and mass (kg/t) material indices, i.e., the volume and max mass of materials processed at the stage per tonne of final product; φmin p,i and φp,i , the lower and upper bounds of the admissible filling degree of the unit; τp,i , the processing time in a vessel; σp,i , the specific capacity of a filter or dryer; hp,i , the fraction of main operations in the total time during which a filter or dryer is occupied by a batch; δp,i , the admissible cake thickness in a filter press; and the durations of loading, unloading, and cleaning of a vessel.

1.2.2

Batch Size and the Compatibility Constraint

Each stage i ∈ Ip with the selected unit x̂i can process a batch of product p only within a certain min ≤ w ≤ w max . For vessels, the bounds of this range are determined by the range of sizes wp,i p p,i 12

admissible filling degrees; for filter presses, by the admissible cake thickness:

min wp,i (x̂i ) =

 min   φp,i x̂i 

up,i vp,i

,

  0,

max wp,i (x̂i ) =

i ∈ A,

(1.7)

i ∈ B ∪ C ∪ D,

 max φp,i x̂i   ,    u v   p,i p,i

δp,i x̂i

,   up,i vp,i      +∞,

i ∈ A ∪ D, (1.8)

i ∈ C, i ∈ B.

For stages of class B, the batch size is not limited directly: the main size of the unit determines not the admissible amount of material but the time needed to process it (see Section 1.2.3). For example, in the CES shown in Fig. 1.1, the product isolation and suspension preparation stages are equipped with vessels (class A), whose main size is the working volume, while the filtration and drying stages (class B) are equipped with units whose main size is the working surface area. This sequence of stages (isolation, filtration, slurrying, drying) forms the final section of the routes of all three products in the test case (stages 122, 134, 138, and 145; Fig. 1.4). A batch of product p must be processed by the units of all stages of the route; therefore, the admissible batch sizes form the intersection of the intervals defined by the individual stages. Consider a fragment of a CES variant comprising stages 1 to i (which corresponds to a node at level i of the search tree, Fig. 1.2) and introduce the quantities min wpmin (Ωe , i) = max wp,l (x̂l ),

wpmax (Ωe , i) =

l∈Ip , l≤i

max min wp,l (x̂l ).

l∈Ip , l≤i

(1.9)

The fragment is operable if the corresponding interval is nonempty for every product. We define the function S(Ωe , i) =

P X 

wpmin (Ωe , i) > wpmax (Ωe , i) , 

(1.10)

p=1

where [ · ] is the indicator function, equal to one if the condition in brackets holds and to zero otherwise. The value S(Ωe , i) equals the number of products for which there is no batch size admissible simultaneously at all stages 1 to i of the fragment. For i = I, the fragment coincides with the complete CES variant; hence the function S(Ωe ) in (1.4) is a special case of the function just introduced: S(Ωe ) = S(Ωe , I). If S(Ωe ) = 0, then constraint (1.4) is satisfied and the CES variant is operable (for every product, there exists a batch size admissible at all stages of its process route), so it makes sense to proceed with checking the processing time constraint (1.5). A nonzero value of S(Ωe ) means that the variant is not operable, and the value itself is the number of products whose batches cannot be processed. Function (1.10) has a property that is fundamental for the branch-and-bound method: when another stage is added to the fragment, the maximum in (1.9) can only increase and the minimum can only decrease, i.e., the interval of admissible batch sizes can only shrink. Consequently, S(Ωe , i) ≤ S(Ωe , i + 1),

i = 1, I − 1,

(1.11)

and if the fragment corresponding to a node of the search tree is not operable, then all CES variants in the subtree of this node are not operable either. It is precisely this property that justifies checking S(W, i) at the internal nodes of the tree in the algorithms considered in the subsequent chapters.

13

If a CES variant is operable, the batch size of each product is set to the largest admissible value: wp = wpmax (Ωe , I), p = 1, P . (1.12) The upper bound is chosen because the number of system cycles needed to produce the required amount of product is inversely proportional to the batch size (see (1.23)), while the processing times in vessels do not depend on it. This makes precise the statement made at the beginning of the chapter that, for an operable CES variant, the operating parameters of the system are determined uniquely. Note that in design problems where the unit sizes are not known in advance, batch sizes are chosen differently, namely, so as to load the equipment of the coupled stages as evenly as possible [58]; in the formulation considered here, where CES variants are composed of standard sizes, no such adjustment is needed. For example, for the solution given in Table 1.2, when the first product is manufactured, the lower bound of the batch size range is determined by stage 138 (slurrying) and the upper bound by stage 115 (diazotization); the range is [2.459, 4.437] t, and according to (1.12) we take w1 = 4.437 t.

1.2.3

Operation Durations

Let dw p,i denote the processing time of a batch of product p in a unit of stage i. For vessels, it is given by the process specification; for filters and dryers, it is determined by the amount of material in the unit and by the unit capacity:   τp,i ,      hp,i gp,i up,i wp    ,   σp,i x̂i 

i ∈ A,

 ,   vp,i σp,i      hp,i µp,i up,i wp    ,

i ∈ C,

dw p,i =  hp,i µp,i δp,i

σp,i F (x̂i )

i ∈ B, (1.13)

i ∈ D,

where gp,i is either the mass material index µp,i or the volumetric material index vp,i (depending on whether the specific capacity σp,i refers to the mass or to the volume of material), and F (x̂i ) is the heat-transfer area of a vacuum dryer with working volume x̂i , taken from the equipment catalogue. For filter presses, the processing time depends neither on the batch size nor on the unit size, since filtration continues until a cake of the given thickness has formed. The cleaning (regeneration) time of filters and dryers is proportional to the duration of the main operations: 1 − hp,i , i ∈ B ∪ C ∪ D; (1.14) dcp,i = dw p,i hp,i for vessels, the durations of loading, unloading, and cleaning are given by the process specification. In the test case (Section 1.3), only the processing times are specified for vessels, and the initial durations of the other operations are set to zero. The actual durations of loading and unloading are then determined when the schedule is computed, as a result of their alignment with the coupled operations of the neighboring stages. The transfer of a portion of material from one stage to another is a pair of simultaneous operations: unloading from the unit of the source stage and loading into the unit of the receiving stage. Both operations of the pair must start and end at the same moments; therefore, their durations are aligned to the larger of the two. If loading at the receiving stage overlaps with processing (for example, a suspension is fed continuously to a filter), unloading from the unit of the source stage continues throughout the filtration.

14

Thus, in the CES of Fig. 1.1, the suspension is unloaded from the unit of the product isolation stage throughout the filtration; hence the duration of this unloading is determined by the filter capacity, i.e., by its working surface area x̂2 and by the batch size. Expressions (1.8), (1.12), and (1.13) show that a CES variant affects the production time in two ways: through the batch sizes (and hence the number of system cycles), which depend on the unit sizes at all stages, and through the processing times at the stages of classes B and D, which depend both on the sizes of the corresponding units and on the batch size.

1.2.4

Schedule of a System Cycle

A system cycle in the production of product p is the time interval from the start of the first operation at the initial stages of the route to the end of the last operation at the final stage. During one system cycle, 1 bp = (1.15) min rp,i i∈Ip

batches of the product are produced, and the units of stage i perform κp,i = rp,i bp cycles of their own. If the ratios rp,i of neighboring stages differ, several loadings or unloadings take place during one unit cycle; for a stage i that receives materials from stage i′ ∈ Pp,i and passes them on to stage i′′ , their numbers are (

Llp,i = max

rp,i′ 1, rp,i

)

(

Lup,i = max

,

rp,i′′ 1, rp,i

)

.

(1.16)

The number of portions transferred from the source stage during a system cycle must equal the number of portions received by the receiving stage: κp,i′ Lup,i′ = κp,i Llp,i ,

i′ ∈ Pp,i .

(1.17)

Let Op denote the set of all operations of all unit cycles at all stages of the route during one system cycle. Each operation o ∈ Op has a type, a duration do , a start time tso , and a finish time tfo = tso + do . Operations are related in two ways: by the process sequence within a unit cycle (for an operation o, let o− denote the immediately preceding operation of the same cycle) and by the material link between stages (for a loading operation o, let o∗ denote the coupled unloading operation at the source stage). The start times of the operations are determined recursively:

tso =

 f  to− , 

ts ∗ ,

o   θ , c

if o− exists, if o is the first in its cycle and o∗ exists, if o is the first in cycle c of an initial stage,

(1.18)

where θc is the start time of cycle c of the units of the initial stage (in the initial schedule, θc = 0). From the data-structure point of view, the schedule of a system cycle is a weighted directed acyclic graph whose vertices are the operations, whose arcs are the relations “o− → o” and “o∗ → o”, and whose vertex weights are the operation durations. The start and finish times of all operations are computed in a single pass over this graph in an order consistent with the direction of the arcs. Figure 1.3 schematically shows the schedule of two consecutive operating cycles of the CES from Fig. 1.1 for rp,i = 1 and νi = 1 (the numbers on the processing bars are the system cycle numbers). At the filtration and drying stages (class B), loading and unloading overlap with processing. Therefore, unloading from the unit of the product isolation stage and loading of the unit of the suspension preparation stage are performed simultaneously with filtration, and unloading of the suspension is performed simultaneously with drying (dashed lines in the figure). As a consequence, the unit of the product isolation stage is occupied longer than its processing operation lasts: its occupation time also includes the filtration time. 15

Figure 1.3: Gantt chart of two CES operating cycles The durations of “filled idle times” are not part of the input data; they are determined by the schedule. A “filled idle time” during loading occurs between two consecutive loadings of a unit if the portions arrive at intervals; its duration equals the gap between the end of one unloading operation and the start of the next one at the source stage and is fully determined by the preceding stages. A “filled idle time” during unloading is more difficult to determine: its duration depends on the operations not only of the immediately following stage but also of all subsequent stages. It is introduced to satisfy the non-overlap condition for cycles: the next cycle of a unit cannot start before the previous one has finished, tsfirst(c′ ) ≥ tflast(c) ,

(1.19)

where first(c) and last(c) are the first and last operations of cycle c, and c′ is the index of the next cycle of the same unit: c′ = c + νi for νi > 1 units operating out of phase and c′ = c + 1 otherwise. The durations of “filled idle times” during unloading are determined by an iterative algorithm that minimizes the duration of the system cycle [58]: 1. The initial schedule (1.18) is constructed with zero durations of the “filled idle times” during unloading. 2. The stages of the route are examined in turn, and condition (1.19) is checked for each cycle c. If it is violated, the overlap ∆ = tflast(c) − tsfirst(c′ ) > 0 is computed. 3. The start of cycle c′ must be delayed by ∆. To this end, moving from the first operation of cycle c′ against the direction of the material flow, one finds the nearest “filled idle time” during unloading at one of the preceding stages after which a portion of material is transferred to cycle c′ , and increases its duration by ∆: the material is held longer in the unit of the source stage. If no such idle time exists and the search has reached the initial stage, the start time θ of the corresponding cycle of the initial stage is increased by ∆. 4. The durations of the “filled idle times” during loading and the start times of the operations are recomputed, after which the check is resumed from step 2. 16

The algorithm terminates when condition (1.19) holds for all cycles at all stages of the route. If rp,i = 1 at all stages of the route, then bp = 1 and κp,i = 1, i.e., during a system cycle each unit performs exactly one cycle with one loading and one unloading; in this case, no “filled idle times” arise and the iterative part of the algorithm is not executed. This is exactly the situation in the test case of Section 1.3. Nevertheless, the software implementation supports the model in full.

1.2.5

Production Time and the Processing Time Constraint

The completed schedule yields the main time characteristics of the CES operating mode in the production of product p. The duration of a system cycle is Tpw = max tfo − min θc . c

o∈Op

(1.20)

Consecutive system cycles overlap: the next cycle starts as soon as the unit that limits the production rate becomes free. The time during which the units of stage i are occupied with the materials of one system cycle is

tp,i =

 κ  p,i  tflast(1) − tsfirst(1) , 

νi > 1, out-of-phase operation,

 tf

otherwise,

νi

(1.21)

s last(κp,i ) − tfirst(1) ,

and the cycle time, i.e., the minimum possible interval between the starts of consecutive system cycles, is determined by the limiting stage: Tpc = max tp,i . i∈Ip

(1.22)

Installing several units operating out of phase at a stage reduces its occupation time by a factor of νi ; in the test case, this applies to stage 138 (slurrying), which is equipped with two units. The number of system cycles required to produce the amount Qp of product p and the production time of this product are ' & Qp Zp = , (1.23) bp wp Tp = Tpw + (Zp − 1) Tpc .

(1.24)

Since the products are manufactured one after another, the production time of the whole product range is the sum of the production times of the individual products, and constraint (1.5) takes the form T (Ωe ) =

P X

Tp ≤ Tmax .

(1.25)

p=1

Unlike the function S, the function T cannot be computed for an incomplete CES variant: the batch sizes (1.12) are determined by the units of all stages of the route, so that until units have been selected for all stages, neither the operation durations (1.13) nor the number of cycles (1.23) are defined. Moreover, T has no monotonicity property similar to (1.11): increasing the unit size at one of the stages may increase the batch size and reduce the number of cycles, but at the same time increase the processing time at the stages of classes B and D and thereby the cycle time. Therefore, constraint (1.5) is checked only for complete CES variants, which are formed at the leaves of the search tree.

17

1.2.6

Optimality Criterion

Let ci,j be the cost of equipping stage i with units of standard size xi,j . The optimality criterion (1.1), i.e., the capital cost of the equipment, is given by Cost(Ωe ) =

I X

ci,ji .

(1.26)

i=1

For a fragment comprising stages 1 to i, the partial sum Cost(Ωe , i) = il=1 cl,jl is nondecreasing in i when the costs are nonnegative; it therefore serves as a lower bound on the cost of all variants in the subtree of the corresponding node. This bound is used to prune branches that cannot contain the optimum (the DefineCost(W, i) procedure in the subsequent chapters). In the test case, the unit costs are taken to be numerically equal to their main sizes; therefore, the value of the criterion for the solution in Table 1.2 equals the sum of the main sizes of the selected units. P

1.2.7

Computational Aspects of the Model

The model determines the computational profile of the problem, which is essential for developing both sequential and parallel algorithms. First, the evaluation of a CES variant consists of three procedures of very different computational cost. Computing the criterion (1.26) and the compatibility function (1.10) amounts to a single pass over the stages of the fragment and can be updated incrementally when moving from a tree node to its child. Computing the function T , in contrast, requires constructing the schedule of a system cycle for each product, i.e., computing the durations and start times of all operations in the sets Op and, in the general case, also the iterative alignment of the cycles. This explains the ratio of the execution times of the procedures S(W, i) and T(W) reported in Chapter 2 (Table 2.1) and justifies the strategy of the sequential algorithm: the expensive check (1.5) is performed only for those complete variants that have passed all inexpensive checks. Second, the structure of the operation graph, i.e., the composition of the sets Op , the operation types, and the relations “o− → o” and “o∗ → o”, is determined only by the process routes and by the parameters rp,i and νi , and does not depend on the CES variant. Only the vertex weights of the graph (the operation durations) and the cycle start times θc depend on the variant. Therefore, it suffices to build the schedule structure once before the search starts, representing it as arrays of indices of the preceding and coupled operations; after that, evaluating each variant amounts to filling the arrays of durations and traversing the graph. This separation of data into an immutable part shared by all variants and a small mutable state specific to a particular variant is of fundamental importance for the GPU implementation: the immutable part can be placed in constant or shared memory and used jointly by all threads, whereas each thread needs its own copy of only the arrays of durations and cycle start times (see Chapters 4 and 5). Third, the cost of computing T is determined by the number of operations in the sets Op , which grows with the number of products, the length of the routes, and the batch size ratios, and, when batches are split and merged, also by the number of cycle alignment iterations, which depends on the input data. The latter makes the computations irregular: the time needed to evaluate different CES variants may vary, which must be taken into account when distributing work among parallel workers. The main notation of the model is summarized in Table 1.1.

1.3

Test Case

The algorithms developed in this book were studied on the equipment selection problem for a CES consisting of 16 processing stages (see Fig. 1.4).

18

Table 1.1: Main notation of the mathematical model of CES operation Symbol

Meaning

p = 1, P Qp Ip , Pp,i A, B, C, D x̂i νi rp,i , up,i vp,i , µp,i max φmin p,i , φp,i min wp ; wp , wpmax S(Ωe , i)

product index production amount of product p, t set of stages in the route of product p; set of source stages for stage i stage classes by the type of the main unit main size of the unit of stage i in variant Ωe number of identical units at stage i batch size ratio; fraction of a batch in one unit volumetric and mass material indices bounds of the admissible filling degree of a unit batch size; bounds of the admissible batch size range, t number of products whose batches cannot be processed by the fragment of stages 1, . . . , i processing and cleaning times, h number of batches per system cycle; number of cycles of the units of stage i per system cycle set of operations of a system cycle; start and finish times of operation o system cycle duration and cycle time, h number of system cycles and production time of product p cost of equipping stage i with units of standard size xi,j

c dw p,i , dp,i bp , κp,i

Op ; tso , tfo Tpw , Tpc Zp , T p ci,j

Within 308 days (7392 h), the CES must produce three dyes in amounts of 2250, 100, and 50 t, respectively. For each stage, a finite set (from 2 to 12) of standard sizes of standard equipment is available for installation. Units of one standard size, chosen from the available set, are installed at each stage. The number of units is given by the input data (two units are installed at stage 138 and one unit at each of the other stages). The problem size in the algorithm tests is set by the number of standard equipment sizes available for each stage. For example, if two standard sizes are available for each stage, the test case shown in Fig. 1.4 has 216 (65 536) possible CES variants. If 5 standard sizes are available, there are 516 (≈ 1011 ) possible CES variants; if 12 standard sizes are available, there are 1216 (≈ 1017 ) possible CES variants. Table 1.2: Solution of the equipment selection problem for the test CES Stage no. 101 103 106 108 110 115 127 125 111 112 118 122 123 134 138 145

Stage

Main size of equipment

Dissolution Diazotization Preparation of the azo component Clarifying filtration Azo coupling Diazotization Dissolution Preparation of the azo component Filtration Dissolution Azo coupling Dye isolation Dissolution Dye filtration Slurrying Dye drying

19

6.3 m3 10 m3 5 m3 5 m2 16 m3 20 m3 0.4 m3 3.2 m3 25 m2 10 m3 32 m3 32 m3 3.2 m3 50 m2 10 m3 20 m3

Figure 1.4: Test chemical-engineering system Since no actual cost data for the equipment were available, the test problems rely on the assumption that the cost of a unit is directly related to its size: in the input data, the cost of each unit is set numerically equal to its main size (in m3 or m2 , depending on the type of equipment). The optimality criterion (1.1) of a CES variant is thus computed as the sum of the main sizes of the units selected for its stages. This simplification is sufficient for studying the algorithms, since the test problems are not intended to model a particular operating plant. An example solution of the equipment selection problem for the test CES is given in Table 1.2. The value of the optimality criterion for this solution is 248.1.

20

Chapter 2

Sequential Branch-and-Bound Algorithm The branch-and-bound method is one of the popular techniques for solving NP-hard optimization problems in various fields, such as combinatorial optimization and artificial intelligence. The method is based on the idea of successively partitioning the set of feasible solutions into subsets [68]. A lower bound on the objective function is then determined for each subset. If the lower bound is not less than the best solution found so far, the subset is discarded, and the process continues until a single solution is obtained. In effect, branch-and-bound searches the solution space using rules that exclude those regions of the search space in which no solution better than the current one can be found. A branchand-bound scheme has three algorithmic components that can be used to fine-tune the behavior of the algorithm: the search strategy, the branching strategy, and the pruning rules. All three components depend on the specific problem being solved and have a considerable effect on the performance of the algorithm [69]. The chapter is organized as follows. Section 2.1 presents the basic version of the sequential branch-and-bound algorithm for the problem. Section 2.2 estimates the execution times of individual parts of the algorithm, which makes it possible to identify the most time-consuming operations. Section 2.3 proposes an optimization of the sequential algorithm based on these estimates.

2.1

Basic Version of the Algorithm

The whole set of variants of a CES consisting of I stages can be represented as a tree with I levels (see Fig. 2.1). Each level i of the tree corresponds to a stage of the CES. Each edge xi,j of the tree corresponds to the selected unit variant j for stage i (one selected element of the set Xi = {xi,j }, i = 1, I, j = Q 1, Ji ). Each internal node ni,k , i = 1, I, k = 1, Ki , Ki = il=1 (Jl ) of the tree corresponds to a fragment of the CES consisting of stages 1 to i. The transition from an internal node ni,k to one of its children (i.e., choosing the main size of a unit and appending it to the initial fragment of the CES variant obtained at the previous levels of the tree) will be called solution construction. A complete path from the root of the tree to one of its leaves represents a CES variant. The possible CES variants are enumerated by a Depth-First Search (DFS) of the tree [70, 71]. For this problem, depth-first search is preferable to Breadth-First Search (BFS), because the depth of the search tree is small, whereas the number of CES variants at each level grows exponentially. With breadth-first search, the amount of main memory needed to store the tree nodes would also grow exponentially.

21

Figure 2.1: Depth-first traversal of the tree of variants At each level of the tree, the main sizes of the units available for this stage of the CES are enumerated. The selected size is appended to the initial fragment of the CES variant obtained at the previous levels of the tree, and the process continues recursively for all operable initial fragments. Thus, at the last level of the tree, a complete operable CES variant is obtained; the recursion stops, and the variant obtained is compared with the current optimal one. If the new variant has a smaller value of the criterion, it replaces the current optimal solution. A complete traversal of the tree of variants (fixing, for each CES stage, a unit size chosen from the set of standard sizes for that stage and forming a complete CES variant), followed by checking constraints (1.4) and (1.5), would require a very large computational effort. To reduce the number of enumerated variants, an algorithm based on the Branch-and-Bound (B&B) method was developed [55, 72]. Branch-and-bound is widely used to solve optimization problems of this kind in various fields (combinatorial optimization, artificial intelligence, etc.) [68], as well as MINLP and MILP problems [73, 74]. The method is based on the idea of successively partitioning the set of feasible solutions into subsets [68] and discarding those subsets that certainly do not contain an optimal solution. To this end, a lower bound on the objective function over a subset is determined. If the lower bound is not less than the best solution found so far, the subset is discarded, and the process continues until a single solution is obtained. This approach is deterministic and yields the globally optimal solution. The pseudocode of the algorithm is shown in Listing 2.1 (most variable names coincide with those used in the problem statement). The traversal of the tree of variants starts with a call of the procedure EnumerateVariants() at the first level (line 1). This procedure is called recursively down to level I (line 4). The current optimal CES variant is stored in the vector Wopt with I elements, each of which holds the selected main size of the unit for the corresponding stage. The working CES variant is stored in a similar vector W. To reduce the enumeration, the algorithm uses the following scheme. At each level i of the tree, the compatibility constraint (1.4) for the joint operation of the CES equipment is checked (procedure S(), line 9); however, it is checked not for the whole CES but only for its initial fragment, which consists of the levels already passed above and the current one (the CES 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

FindSolution () { E n u m e r a t e V a r i a n t s (1) ; } /* r e c u r s i v e depth - first t r a v e r s a l of the tree */ EnumerateVariants (i) { if ( i <= I ) { for ( j = 1; j <= J [ i ]; j ++) { /* build the initial f r a g m e n t of the CES ( stages 1 to i ) */ W [ i ] = X [ i ][ j ]; /* check the c o m p a t i b i l i t y c o n s t r a i n t (1.4) */ if ( S (W , i ) == 0 && /* cost of the initial f r a g m e n t of the CES ( stages 1 to i ) */ DefineCost (W , i ) < minCost ) { /* r e c u r s i v e search */ E n u m e r a t e V a r i a n t s ( i + 1) ; } } } else { /* all tree levels have been passed ( i == I +1) */ /* check the p r o c e s s i n g time c o n s t r a i n t (1.5) */ if ( T ( W ) <= Tmax ) { /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost (W , I ) < minCost ) { /* store the new optimal s o l u t i o n */ Wopt = W ; minCost = DefineCost (W , I ) ; } } } }

Listing 2.1: Sequential version of the algorithm

fragment consisting of stages 1 to i). If the fragment is already inoperable, there is no point in enumerating units at the lower levels. In addition, the cost (1.1) of the initial fragment of the CES variant (procedure DefineCost(), line 11) is compared with the cost minCost of the current optimal CES variant. Thus, the algorithm descends to the next recursion level (a call of EnumerateVariants() at level i + 1, line 13) only if the initial fragment of the CES consisting of stages 1 to i is, first, operable and, second, not more expensive than the best CES variant found so far. Otherwise, the algorithm backtracks to the previous level of the tree. As a result, the number of enumerated variants is reduced, and the higher the level of the tree at which such a mismatch is detected, the greater the reduction. If all levels of the tree have been passed (i==I+1, line 14), the processing time constraint (1.5) is checked (procedure T(), line 16). If constraint (1.5) is satisfied, a new operable CES variant has been found, and its cost (the result of DefineCost()) is compared with the cost minCost of the current best CES variant (line 18). If a better solution has been found, it replaces the current optimal solution Wopt, and its cost is stored in minCost (line 20). The initial value of minCost is the sum of the maximum unit costs over all stages (the cost of a CES variant consisting of the most expensive units is the largest possible).

2.2

Execution Times of Individual Parts of the Algorithm

For the computational experiments, the algorithm presented above (see Listing 2.1) was implemented in C++. The experiments were carried out on a computer with two Intel Westmere X5650 processors (six cores at 2.6 GHz) and 4 GB of RAM running Red Hat Linux. The program was compiled with Intel C++ Compiler v. 11.1. The test case was the equipment selection problem for the CES described in Section 1.3, with 5 standard equipment sizes available for each processing stage. To identify the main bottlenecks of the algorithm, the execution times of its individual elements were measured. The results are given in Table 2.1. The times in Table 2.1 are negligible for a single evaluation. However, since the search repeats these operations billions of times, the times accumulate and may add up to tens or hundreds of hours. In particular, for the example considered (16 stages with 5 standard sizes available for each), the total run time of the program was 27 h 11 min. The results show that the most expensive operation is the procedure T(W), which checks the

23

Table 2.1: Execution times of individual elements of the algorithm Element of the algorithm Recursive call of the procedure EnumerateVariants(i) S(W, i) T(W) DefineCost(W, i)

Execution time, µs 0.1 4.0 417.0 0.7

processing time constraint. To reduce the total run time of the algorithm, the number of calls of this procedure must be minimized; to this end, the algorithm was optimized as described below.

2.3

Optimization of the Sequential Algorithm

A call of the procedure DefineCost(W), which evaluates the objective function, is a relatively cheap operation (see Table 2.1). If the optimality criterion is checked immediately after an operable CES variant has been obtained, the processing time constraint T(W) need not be checked for solutions that are certainly not optimal. Thus, although the number of calls of DefineCost(W) increases, the number of calls of T(W) decreases noticeably, which ultimately reduces the total run time of the algorithm. The changes made to the algorithm are shown in Listing 2.2. 13 14 15 16 17 18 19 20 21

... else { /* all tree levels have been passed ( i == I + 1) */ /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost (W , I ) < minCost ) { /* check the p r o c e s s i n g time c o n s t r a i n t (1.5) */ if ( T ( W ) <= Tmax ) { /* store the new optimal s o l u t i o n */ Wopt = W ; minCost = DefineCost (W , I ) ; } } } ...

Listing 2.2: Modification of the sequential algorithm The optimality criterion is checked (line 16) before the processing time constraint (line 18). As a result, the run time was reduced approximately by half (13 h 52 min instead of 27 h 11 min).

24

Chapter 3

Parallel Branch-and-Bound for Multi-Core and Cluster Systems Despite the optimization of the algorithm described in Section 2.3, the run time remains considerable. To reduce it further, it is therefore worthwhile to develop and study a parallel version of the algorithm. The tree-like organization of the search space in branch-and-bound is well suited for this, since different branches of the tree can be processed simultaneously [29, 75]. The first results on parallelizing the algorithm for the CES equipment selection problem were reported in [76]. The chapter is organized as follows. Section 3.1 reviews methods for parallelizing search tree traversal: static and dynamic tree partitioning, dynamic load balancing methods and the strategies they use, speedup anomalies, the assigned subtrees and allocated subtrees methods, and approaches to implementing branch-and-bound on GPUs; at the end of the section, the implementations developed in this book are placed in this context. Section 3.2 presents a parallel implementation of the algorithm for shared-memory systems based on the OpenMP standard, in which the nodes of the upper levels of the search tree are treated as independent tasks, and studies the effect of the granularity parameter on the run time. Section 3.3 presents an implementation for distributed-memory systems based on the MPI library and the master–worker paradigm [16, 77]. For both implementations, the speedup and the processor efficiency are studied.

3.1

Methods for Parallelizing Search Tree Traversal

3.1.1

Features of Parallel Search Tree Traversal

Many discrete optimization problems arising in planning, scheduling, computer-aided design, and robotics are formulated as a search in a graph or tree: one has to find a minimum-cost path from the initial node to one of the goal nodes [29, 78]. In general form, the problem is stated as follows: given a finite set X of feasible solutions and an objective function f (x) defined on it, find an element xopt such that f (xopt ) = min{f (x) | x ∈ X}. Constraint satisfaction problems, in which it suffices to find at least one element of a subset S ⊂ X, fit the same scheme if we set f (x) = 0 for x ∈ S and f (x) = 1 for the other elements of X. In practically important problems, the cardinality of X is so large that exhaustive enumeration of its elements is impossible. Branch-and-bound, dynamic programming, and heuristic search exploit the structure of the tree (graph) of variants in order to exclude entire branches from consideration and find a solution without explicitly enumerating all elements of X [29, 69]. The optimal CES equipment selection problem considered in Chapter 1 belongs to exactly this class: each level of the tree corresponds to a CES stage, each leaf to a complete CES variant, and checking the constraints and the cost of initial fragments makes it possible to prune subtrees 25

that certainly do not contain the optimal solution (see Section 2.1). A search tree naturally decomposes into independent subtrees that can be traversed simultaneously. This makes tree traversal methods a convenient target for parallelization [75, 79]. In branch-and-bound, however, the tree is not given in advance: it is built during the computation, and its shape depends on the input data and on how early a good solution allowing branches to be pruned is found. It is therefore unknown in advance how many nodes will be visited in each subtree, and the central problem of parallel traversal is to distribute the work among processors so as to minimize the total traversal time (or, equivalently, to maximize the efficiency of using the computing system) [80, 81]. Efficiency losses in parallel tree traversal consist of three components: • idle time of processors that have received smaller subtrees than the others (load imbalance); • communication overhead: transferring subtrees, work requests, and solutions found; • redundant work, i.e., nodes that are visited by the parallel algorithm but would not be visited by the sequential one (for details, see Section 3.1.5). Different parallel traversal methods trade off these components in different ways, and the choice of method depends both on the properties of the tree of the problem at hand and on the architecture of the computing system.

3.1.2

Static and Dynamic Tree Partitioning

According to how the tree is split into subtrees and how these are distributed among processes, all methods of parallel depth-first tree traversal can be divided into two groups: those with static and those with dynamic partitioning [78, 81]. With static partitioning, the tree is split into subtrees once, before the traversal begins, and each subtree is assigned to a specific process. For example, each process can uniquely compute the roots of “its” subtrees from its rank and need not interact with the other processes afterwards. This approach is simple to implement and requires practically no communication; however, it is efficient only when the subtrees are known to be of similar size, i.e., for trees with a regular structure known in advance. Similarly, static decomposition of the computational domain is efficient in problems with a uniformly distributed computational load, for example, in solving the unsteady heat equation with an explicit finite-difference scheme [82]. For branch-and-bound, this condition usually does not hold: pruning makes the tree irregular, so that some subtrees are traversed in fractions of a second while others take hours. With dynamic partitioning, the distribution of work changes during the traversal. A process that has finished traversing its part of the tree becomes idle and receives part of the notyet-traversed work of one of the busy processes. The process that gives away work will be called the donor, and the process that receives it, the receiver. Dynamic partitioning keeps the processors evenly loaded for trees of arbitrary shape at the price of additional communication and task-scheduling overhead. This is why dynamic schemes are predominantly used for branchand-bound [79, 83].

3.1.3

Classification of Dynamic Load Balancing Methods

Dynamic load balancing methods differ in who initiates the splitting of the tree into subtrees and who initiates the transfer of subtrees to other processes. Three groups of methods are distinguished [80, 81] (Table 3.1). In the methods of the first group, the donor sends the nodes it generates to other processes regardless of their state. In the simplest variant, each new node is sent to a randomly chosen process [84]; on average, this balances the load, but it generates communication at every branching step and overloads the network. Moreover, the number of nodes accumulated by an individual 26

Table 3.1: Classification of dynamic search tree partitioning methods Initiator of splitting and transfer

Typical scheme

Advantages and drawbacks

Splitting and transfer are initiated by the donor

Generated nodes are sent to randomly chosen or neighboring processes

No requests needed; large communication volume, unbounded growth of the number of nodes held by a process

Splitting is initiated by the donor, transfer by an idle process

Master processes generate a pool of subtree roots in advance; workers obtain them on request (master–worker)

Simple control, low communication volume; possible imbalance due to different subtree sizes

Splitting and transfer are initiated by an idle process

An idle process requests work; a busy process hands over part of its subtree

Best balancing for irregular trees; costs of requests, rejections, and work splitting

process is unbounded, which may lead to memory shortage. The communication volume is reduced if nodes are sent only to neighboring processors or if nodes below a certain cutoff level may not be transferred and are processed by the process that generated them [85]. In the methods of the second group, work is generated in advance but handed out to idle processes only on request. The processes are divided into masters, which generate subproblems, and workers, which solve them. A typical representative of this group is the master–worker paradigm, in which a single master process traverses the upper part of the tree and forms a pool of subtree roots, and the worker processes receive them one at a time as they become free. To reduce the load on the master process, schemes with several sources of subproblems are used, each serving its own group of workers; when a source runs out of subproblems, its workers switch to another source. A similar idea underlies the fixed-packet DFS method [86]: all processes first traverse the tree synchronously down to a certain level, the resulting nodes are divided into as many classes as there are processes, and then processes that have become idle request some of the root nodes from the classes of other processes. The third group, which is the most widely used in practice, comprises methods in which both the splitting and the transfer of work are performed at the request of an idle process [80, 85, 87, 88]. Each process traverses its assigned subtree depth-first and periodically checks whether requests from other processes have arrived. Having finished its traversal, a process selects a potential donor and sends it a request. The donor either rejects the request if its remaining work is too small, or splits its work into two parts and sends one of them to the receiver. In modern terminology, this approach is called work stealing. For multithreaded shared-memory computations, it has been proved to provide near-optimal load balance with moderate overhead [89].

3.1.4

Algorithmic Aspects of Dynamic Load Balancing

The methods of the third group differ in the strategies used at each step of work exchange. Donor selection. With random polling, an idle process chooses a donor at random and repeats the choice if it is rejected [80, 88]. This strategy scales well, but when little work remains, it leads to many rejections. With a global round robin, all processes access a common pointer, which is advanced to the next process after each request. This considerably reduces the number of requests for trees of various structures; however, with many processes, contention for the shared pointer arises and it becomes a bottleneck [80]. Finally, with targeted selection, a master process that receives information from the workers about the traversed parts of their subtrees directs the request to the most heavily loaded process [81]. Work splitting strategy. The donor may give the receiver a single subtree by removing one 27

not-yet-visited node from its stack (node splitting), or split its stack into two parts so that each part contains approximately the same number of nodes of different levels (stack splitting) [80, 85]. The second option transfers more work per exchange and therefore reduces the number of requests, but requires more data to be sent. To transfer as large a subtree as possible, the unvisited node closest to the root is chosen. Initial work distribution. The traversal may start with a single process to which the whole tree is assigned; the other processes then obtain work through the request mechanism [85, 88]. This is simple, but at the beginning most processes are idle waiting for work, and the number of interprocess exchanges rises sharply. An alternative is a preliminary traversal of the upper part of the tree (usually breadth-first) until enough root nodes have been obtained, which are then distributed among the processes, or having each process compute the root of its first subtree from its own rank. Splitting threshold. Splitting very small subtrees is not worthwhile, because the communication overhead would exceed the gain from parallel traversal. Therefore, the donor rejects a work request if, for example, the current traversal depth is close to the maximum or the subproblem contains too few nodes [81]. The specific threshold depends on the problem size, the cost of processing a node, and the characteristics of the communication network. Using these strategies, parallel algorithms have been developed for the travelling salesman problem, the knapsack problem, the Boolean satisfiability problem, puzzles, and a number of other combinatorial problems. For most of them, near-linear speedup has been achieved on multiprocessor systems [78, 79].

3.1.5

Speedup Anomalies and the Search Overhead Factor

A specific feature of parallel search is that the parallel and the sequential algorithms usually visit different sets of nodes. If any solution is to be found (as in constraint satisfaction problems), the search stops as soon as one of the processes finds a feasible solution, and the parallel algorithm may visit either fewer or more nodes of the search tree than the sequential one. In branch-andbound, the same effect arises because the incumbent (the best value found so far), which is used to prune branches, is updated at different moments. If a good solution is found by one of the processes earlier than the sequential algorithm would find it, the other processes prune more branches; if later, part of the work turns out to be redundant. To quantify this effect, the search overhead factor is introduced: R=

Wp , Ws

(3.1)

where Ws is the number of nodes visited by the sequential algorithm and Wp is the total number of nodes visited by all processes of the parallel algorithm [29, 90]. The speedup on s processors is then bounded from above by s/R. For R > 1, the parallel algorithm performs redundant work, and the speedup is necessarily sublinear. For R < 1, superlinear speedup is possible: the parallel algorithm visits fewer nodes than the sequential one. Such situations, known as speedup anomalies, have been studied in detail for parallel branch-and-bound [91]. It has been shown that increasing the number of processors may either accelerate or slow down the search, and conditions have been formulated under which a slowdown is impossible. For parallel depth-first search with backtracking, it was shown in [90] on abstract models of the search space that if solutions are distributed identically over the search regions of different processes, then the amounts of work of the sequential and parallel algorithms coincide on average. If the solution distributions differ, the sequential algorithm visits on average more nodes than the parallel one. A systematic superlinear speedup rather indicates that the chosen sequential traversal order is not the best one for the given problem.

28

3.1.6

The Assigned Subtrees and Allocated Subtrees Methods

For traversing the trees of combinatorial problems on distributed-memory multiprocessor systems, two methods were proposed in [81, 92]; they represent the second and the third group of the above classification, respectively. In both methods, the processes are divided into one master process and several workers. The assigned subtrees method focuses on reducing interprocess communication. The master process traverses the tree breadth-first from the root until the number of generated but not yet processed nodes (subtree roots) reaches a given value m, which considerably exceeds the number of worker processes. The root nodes are then handed out to the workers one at a time as they finish traversing the subtrees received earlier. The workers interact with the master process only when receiving a new task and sending back a result, so the control overhead is small. However, the method does not guarantee a uniform load: individual subtrees may be many times larger than the others, and at the end of the computation some processes will be idle. Increasing m reduces the size of individual subtrees and smooths out the imbalance, but increases the number of exchanges with the master process. Thus, the assigned subtrees method is advisable when the partitioning yields subtrees of similar size or when there are enough of them for the differences in size to average out. The allocated subtrees method emphasizes uniform load. The master process generates only as many root nodes as there are workers and hands them out one per worker. When a worker finishes its traversal, the master process chooses, among the busy processes, the one with the largest untraversed part of its subtree (the workers periodically report the progress of their traversal to the master) and instructs it to split off a subtree for the idle process. A subtree is split off only if a splitting condition holds, for example, if the index of the last fully traversed level is much smaller than the total number of tree levels. The allocated subtrees method ensures a uniform load, but requires additional effort for collecting progress information and for splitting off subtrees.

3.1.7

Parallelizing Branch-and-Bound on GPUs

The methods discussed above were developed for multiprocessor (multi-core) systems in which each processor traverses its own subtree independently. The GPU architecture imposes additional constraints on the parallelization of branch-and-bound. Thousands of GPU threads are grouped into warps that execute the same instruction; the branching inevitable in the traversal of an irregular tree causes different code paths within a warp to be executed serially (divergence). The amount of fast memory per thread is small, and dynamic memory allocation and data exchange between threads of different blocks are difficult. Finally, host–device transfers are relatively expensive, so frequent movement of work between the CPU and the GPU is undesirable. Existing approaches to implementing branch-and-bound on GPUs can be roughly divided into three groups. In the first group, the search tree is managed on the CPU, while the most computationally intensive part, the simultaneous processing (branching and bounding) of a large number of nodes from the pool of active nodes, is offloaded to the GPU [33, 93]. In the second group, the tree is partitioned in advance into many subtrees, each of which is traversed entirely by a separate GPU thread; the implementations presented in Chapter 4 also belong to this group [94, 95]. In the third group, both the traversal and the dynamic load balancing are implemented on the GPU; for example, the Integer-Vector-Matrix data structure makes it possible to represent subtrees compactly as intervals and to redistribute work among threads directly on the accelerator [96]. A comparison of branch-and-bound implementations on multi-core CPUs and on massively parallel accelerators for large problems is given in [97].

29

3.1.8

Methods Used in This Book

The above classification makes it possible to position the parallel branch-and-bound implementations presented in this chapter and in Chapter 4. The implementation for shared-memory systems (Section 3.2) uses the task mechanism of the OpenMP standard: each node of the upper levels of the tree, down to level G, is turned into an independent task, and the OpenMP runtime distributes the tasks among the threads. This is a dynamic scheme in which the task pool serves as the work pool, and the granularity parameter G determines the number and size of the generated tasks. The implementation for distributed-memory systems (Section 3.3) is based on the master– worker paradigm and the assigned subtrees method: the master process traverses the tree down to level G and hands out initial fragments of CES variants to the workers on their requests. In terms of the assigned subtrees method, the number of assigned subtrees m equals the number of Q nodes at level G, i.e., m = KG = G l=1 Jl ; as in the general case, the choice of G determines the trade-off between load uniformity and the amount of communication with the master process. The GPU implementation (Chapter 4) uses static partitioning: the host traverses the tree down to level G, and each of the KG resulting subtrees is assigned to a separate GPU thread without subsequent redistribution of work. This approach minimizes communication and simplifies the kernel code, but makes the run time dependent on the largest subtree, which, together with thread divergence, limits the achievable speedup (see Section 4.4).

3.2

Shared-Memory Systems (OpenMP)

3.2.1

Program Design

In shared-memory systems, several processors have equal access to memory. Such systems are also called Symmetric Multiprocessing (SMP) systems. In the proposed version of the algorithm, all tree nodes Ni = {ni,1 , ni,2 , . . . , ni,k }, i = 1, I, k = Q 1, Ki , Ki = il=1 Jl at each level i are treated as independent tasks that can be executed in P parallel. The total number of such tasks, Ntasks = G i=1 Ki , 1 ≤ G ≤ I, may be very large. The granularity parameter G of the algorithm limits the degree of parallelism to a certain level of the tree: subtrees below this level are not split into parallel tasks but are processed sequentially [54]. Tasks are an alternative to traditional threads; they provide a high level of abstraction and allow the programmer to express the logical parallelism of an application without managing parallel execution directly [98, 99]. Tasks provide better load balancing and efficient use of the available computing resources. The pseudocode of the proposed parallel algorithm is shown in Listing 3.1. The parallel implementation is a development of the sequential program considered above (see Chapter 2, Listings 2.1 and 2.2); directives of the Open Multi-Processing (OpenMP) standard [100] were added to its source code. Starting with version 3.0, OpenMP supports the task directive, which makes it possible to create separate independent tasks. The initial call of the recursive procedure EnumerateVariants() (line 4 of Listing 3.1) is made inside a #pragma omp single directive nested in a #pragma omp parallel directive, so that only one of the available parallel threads performs the recursive tree traversal. Subsequent recursive calls of EnumerateVariants() (line 16) are made inside the #pragma omp task directive. Thus, while recursively traversing the tree, the first thread of the program creates new independent tasks at each level i, 1 ≤ i ≤ G; these tasks are automatically placed into the task pool, from which any idle thread can take them for execution. Note that, when these tasks are processed in parallel, each task must be given its own copy of the array W, the current initial fragment of the CES.

30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

FindSolution () { # pragma omp p a r a l l e l { # pragma omp single { E n u m e r a t e V a r i a n t s (1) ; } } } EnumerateVariants (i) { if ( i <= I ) { for ( j =1; j <= J [ i ]; j ++) { /* build the initial f r a g m e n t of the CES ( stages 1 to i ) */ W [ i ] = X [ i ][ j ]; /* check the o p e r a b i l i t y and the cost of the initial f r a g m e n t of the CES */ if ( S (W , i ) == 0 && DefineCost (W , i ) < minCost ) { if ( i <= G ) { /* current search depth does not exceed G */ /* create a p a r a l l e l task */ # pragma omp task E n u m e r a t e V a r i a n t s ( i + 1) ; } else { /* current level is deeper than G */ /* s e q u e n t i a l r e c u r s i v e tree t r a v e r s a l */ E n u m e r a t e V a r i a n t s ( i + 1) ; } } } } else { /* all tree levels have been passed ( i == I + 1) */ ... } }

Listing 3.1: Parallel version of the algorithm for shared-memory systems

3.2.2

Results of Computational Experiments

For the computational experiments, a parallel version of the program was developed in C++ using OpenMP 3.0. The test case was the equipment selection problem for the CES described in Section 1.3, with 5 standard equipment sizes available for each processing stage. The experiments were carried out on a computing system with four eight-core Intel Xeon E7550 processors at 2 GHz (32 cores in one machine) and 128 GB of RAM running Red Hat Linux. The program was compiled with Intel C++ Compiler 11.1. Figure 3.1 shows how the run time of the parallel program on different numbers of processor cores depends on the granularity parameter G. For values of G greater than 10, the number of tasks became so large that there was not enough memory to store them (the program terminated abnormally). The speedup and the efficiency of using the cores were also studied. In general, the speedup Sp is defined as the ratio of the time Ts needed to execute the sequential program on one core to the time Tp needed to execute the parallel program on a p-core machine: Sp = Ts /Tp . The case Sp = p corresponds to linear (ideal) speedup. The efficiency Ep is defined as the ratio of the speedup Sp to the number of cores p: Ep = Sp /p = Ts /(p · Tp ). The results are shown in Fig. 3.2 (the horizontal axis is logarithmic). The granularity parameter was G = 10. It is recommended to choose the granularity parameter G so that the number of tasks considerably exceeds the number of physical cores of the computing device and is at least not smaller than the number of cores used. The number of tasks is determined by the value of G, for which we propose the empirical estimate G ≥ logmax{Ji } Np , i = 1, I, where max{Ji } is the maximum cardinality of the set of unit sizes available at a stage and Np is the number of cores used. However, with a very large number of tasks (in our case, for G > 10), memory may run out. It is difficult to determine the required amount of memory exactly. Nevertheless, the results in Fig. 3.1 show that, starting from a certain threshold value of G (in our example, G = 4), the run time of the program remains practically unchanged despite the increasing number of tasks. We therefore suggest first performing computations for a small problem, finding the threshold value of G that yields an acceptable solution time, and then carrying out the computations for larger problems. The study of speedup and efficiency showed that the parallel program achieves near-linear 31

25000

Run time, s

20000

4 cores 8 cores 16 cores 32 cores

15000 10000 5000 0

1

2

4

6

Granularity parameter G

8

10

Figure 3.1: Run time of the parallel program versus the granularity parameter G speedup and fairly high scalability. The efficiency of using the cores is 0.9–0.95.

3.3

Distributed-Memory Systems (MPI)

In distributed-memory systems, each compute node exchanges data with the other nodes only over communication links, by message passing. The parallel implementation of the algorithm for distributed-memory systems is based on the master–worker paradigm (also known as the master–slave paradigm) [101, 102] and on the assigned subtrees method [92]. The main features of the master–worker paradigm are the following: • the master process (master) generates tasks and distributes them evenly among the worker processes (workers); • after solving the task it has received, each worker returns the results to the master and requests a new task; • the workers do not exchange data with one another, only with the master; • after receiving the results of a task, the master dynamically generates a new task for the worker that has become free or, if no tasks remain, sends it a termination signal; • the algorithm requires at least two processes (one master and one worker). The pseudocode of the master and worker processes is shown in Listings 3.2 and 3.3, respectively.

3.3.1

Master Process

Using the recursive procedure MasterEnumerateVariants(), the master performs a sequential depth-first traversal of the tree only down to a certain level G (the granularity), 1 ≤ G ≤ I, forming the initial fragments of CES variants W [i], i = 1, G (lines 23–29). 32

35

1

30

0.8

25

Linear speedup Measured speedup Efficiency

15

Efficiency

Speedup

20

0.6

0.4

10

0.2

5 0

1

2

4

8

Number of processors

16

32

0

Figure 3.2: Speedup and efficiency of using the cores (G = 10) The granularity parameter G is used to vary the degree of parallelism in the same way as in the version of the algorithm for shared-memory systems (see Section 3.2): subtrees below this level are not split into parallel tasks but are processed sequentially. Its value must be such that the number of generated subproblems is at least not smaller than the number of processing elements used. At the lowest level of recursion, the master waits for messages from the workers (line 33); these messages can be of two types: solution found (SOLUTION_FOUND) or work request (REQUEST_WORK). If the master receives a new complete CES variant from a worker (line 34), its cost is compared with the cost of the current global optimal CES variant (line 36). If a new optimal CES variant has been found, it replaces the current optimal solution Wopt, and its cost is stored in minCost (line 38). If the master receives a request for a new task from an idle worker (line 39), it sends this worker a do work message (DO_WORK) containing the initial fragment of a CES variant and the cost of the current optimal CES variant (line 41). The variable minCost is initialized with the sum of the maximum unit costs over all stages.

3.3.2

Worker Process

A worker requests a task from the master (line 4) and waits for a reply. The reply can be of two types: do work (DO_WORK) or quit (QUIT). If a message with a new task has been received from the master (line 7), containing an initial fragment of a CES variant and the cost of the current optimal CES variant, the recursive procedure WorkerEnumerateVariants() is called (line 10). Using this procedure, the worker performs a sequential depth-first traversal of the subtrees assigned to it, from level G + 1 down to the lowest level I, completing the initial fragment received from the master to a full CES variant W [i], i = G + 1, I (lines 20–27). The transition to the next recursion level is made in the same way as in the sequential version of the algorithm: only if the initial fragment of the CES consisting of stages 1 to i is, first, operable and, second, not more expensive than the best CES variant found so far (line 25). Otherwise, the algorithm backtracks to the previous level of the tree. When all levels of the 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

Master () { /* number of workers ( one of the p r o c e s s e s is the master ) */ num_workers = NUM_PRO CESSORS - 1; /* start the tree t r a v e r s a l */ M a s t e r E n u m e r a t e V a r i a n t s (1) ; /* while some workers are still running */ while ( num_workers > 0) { msg = R e c e i v e W o r k e r M e s s a g e () ; /* a worker has sent a new s o l u t i o n */ if ( msg . type == SOLUT ION_FOUN D ) { /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost ( msg .W , I ) < minCost ) { /* store the new optimal s o l u t i o n */ Wopt = msg . W ; minCost = DefineCost ( msg .W , I ) ; } } /* a worker r e q u e s t s a new task */ elseif ( msg . type == REQUEST_WORK ) { /* stop the worker */ S e n d W o r k e r M e s s a g e ( msg . workerID , QUIT ) ; num_workers - -; } } } MasterEnumerateVariants (i) { if ( i <= G ) { /* current depth does not exceed G */ for ( j =1; j <= J [ i ]; j ++) { /* build the initial f r a g m e n t of the CES ( stages 1 to i ) */ W [ i ] = X [ i ][ j ]; /* check the o p e r a b i l i t y and the cost of the initial f r a g m e n t */ if ( S (W , i ) == 0 && DefineCost (W , i ) < minCost ) { M a s t e r E n u m e r a t e V a r i a n t s ( i + 1) ; } } } /* r e c u r s i v e search */ else { /* an initial f r a g m e n t of G stages has been built ( i == G + 1) */ while ( true ) { /* receive m e s s a g e s from workers */ msg = R e c e i v e W o r k e r M e s s a g e () ; if ( msg . type == SOLUT ION_FOUN D ) { /* a s o l u t i o n r e c e i v e d */ /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost ( msg .W , I ) < minCost ) { /* store the new optimal s o l u t i o n */ Wopt = msg . W ; minCost = DefineCost ( msg .W , I ) ; } } elseif ( msg . type == REQUEST_WORK ) { /* work request */ /* send the initial f r a g m e n t of the CES and the current c r i t e r i o n value */ S e n d W o r k e r M e s s a g e ( msg . workerID , DO_WORK , W , minCost ) ; break ; } } } }

Listing 3.2: Pseudocode of the master process

tree have been passed (i==I+1, line 28), the cost of the variant is compared with the cost of the current local best CES variant (line 30). If a better solution has been found, the processing time constraint (1.5) is checked (procedure T(), line 32). If constraint (1.5) is satisfied, a new optimal CES variant has been found; it replaces the current local optimal solution Wopt, and its cost is stored in minCost (line 34). After the recursive procedure has finished, the worker sends the local optimal CES variant to the master (line 13). The worker then becomes free and sends the master a request for a new task. On receiving a termination message from the master (line 15), the worker terminates (line 16).

3.3.3

Results of Computational Experiments

For the computational experiments, a program was developed in C++ using the Message Passing Interface (MPI) library [103]. The test case was the equipment selection problem for the CES described in Section 1.3, with 5 standard equipment sizes available for each processing stage. The program was tested on a computing cluster of 198 nodes, each equipped with two six-core processors (Intel Westmere X5650, 2.6 GHz, 2 or 4 GB of RAM) running Red Hat Linux. The cluster nodes are connected by a high-speed, low-latency switched InfiniBand network. The 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Worker () { while ( true ) { /* request a new task */ S e n d M a s t e r M e s s a g e ( workerID , REQUEST_WORK ) ; msg = R e c e i v e M a s t e r M e s s a g e () ; /* a new task has been r e c e i v e d */ if ( msg . type == DO_WORK ) { minCost = msg . minCost ; f o u n d N e w S o l u ti o n = false ; W o r k e r E n u m e r a t e V a r i a n t s ( G + 1) ; if ( f o u n dN e w S o l u t i o n ) { /* a new s o l u t i o n has been found */ /* send the s o l u t i o n to the master */ S e n d M a s t e r M e s s a g e ( workerID , SOLUTION_FOUND , Wopt ) ; } } /* a t e r m i n a t i o n message has been r e c e i v e d */ elseif ( msg . type == QUIT ) { break ; } } } WorkerEnumerateVariants (i) { if ( i <= I ) { for ( j =1; j <= J [ i ]; j ++) { /* build a c o m p l e t e CES variant */ W [ i ] = X [ i ][ j ]; /* check the o p e r a b i l i t y and the cost of the CES variant */ if ( S (W , i ) == 0 && DefineCost (W , i ) < minCost ) { /* r e c u r s i v e search at the next level */ W o r k e r E n u m e r a t e V a r i a n t s ( i + 1) ; } } } else { /* all tree levels have been passed ( i == I + 1) */ /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost (W , I ) < minCost ) { /* check the p r o c e s s i n g time c o n s t r a i n t (1.5) */ if ( T ( W ) <= Tmax ) { /* store the new optimal s o l u t i o n */ Wopt = W ; minCost = DefineCost ( Wopt , I ) ; f o u n d N e w S o l ut i o n = true ; } } } }

Listing 3.3: Pseudocode of the worker process

program was compiled with Intel C++ Compiler 11.1. The parallel program requires at least two processors (one for the master and one for a worker). Figure 3.3 shows how the run time of the parallel program on different numbers of processors depends on the granularity parameter G. The speedup and the processor efficiency were also studied. Figure 3.4 (the horizontal axis is logarithmic) shows the speedup of the program and the processor efficiency. The study of speedup and efficiency showed that the parallel program achieves near-linear speedup and fairly high scalability. However, Fig. 3.4 shows that the processor efficiency and the speedup decrease when the number of processors exceeds 768 (64 cluster nodes). This happens because the master process becomes a bottleneck of the program when it serves a large number of workers, which reduces the speedup and efficiency. This problem can be addressed by the hierarchical master–worker paradigm, in which several submasters are added to the main master process [104, 105].

35

25000

Run time, s

20000

4 nodes 8 nodes 16 nodes 32 nodes 64 nodes

15000 10000 5000 0

1

2

4

6

8

10

Granularity parameter G

12

14

16

Figure 3.3: Run time of the program versus the granularity parameter G

1600

1

1400

0.8

1000

Linear speedup Measured speedup Efficiency

800 600

Efficiency

Speedup

1200

0.6

0.4

400

0.2

200 0

1

2

6

12

24

48

96

Number of processors

192

384

768

0 1536

Figure 3.4: Speedup and efficiency of the program versus the number of processors (G = 10)

36

Chapter 4

Branch-and-Bound on Graphics Processing Units Chapter 3 considered the parallelization of branch-and-bound for shared-memory and distributedmemory systems. This chapter deals with its implementation for heterogeneous computing systems equipped with graphics processing units (GPUs). Modern high-performance computing systems have a heterogeneous architecture that includes a general-purpose multi-core Central Processing Unit (CPU) and a set of various accelerators [106]. The performance of heterogeneous systems is increased not only by increasing the number of processors and cores, but also by incorporating (either integrated or discrete) specialized accelerators for parallel computing, such as Graphics Processing Units (GPUs), Digital Signal Processors (DSPs), and Field-Programmable Gate Arrays (FPGAs) [36]. Such platforms can considerably improve performance in applications such as image processing, medical visualization, and data acquisition, while keeping power consumption fairly low and at a relatively small additional cost. Today, GPUs are among the most powerful computing devices available at an affordable price [107]. Architecturally, a CPU consists of a few cores optimized for sequential processing, whereas a GPU consists of thousands of smaller and more energy-efficient cores designed to handle many tasks simultaneously. Originally, graphics accelerators were devices intended for displaying two- or three-dimensional graphics on a computer screen. As they evolved, they acquired support for so-called shaders, small data-processing programs executed on the graphics accelerator itself. The originally fixed algorithm for computing lighting and transforming vertex coordinates was replaced by a user-defined one. Later, it became possible to write programs computing the color of a pixel on the screen. For this reason, programs for GPUs came to be called vertex and pixel shaders. Shaders enabled developers to create, for example, their own special effects, which considerably increased the complexity and realism of computer graphics and, accordingly, the attractiveness of the resulting images. Shaders had to be programmed in the GPU assembly language, their length did not exceed 20 instructions, branch instructions were not supported, and computations were performed in fixed-point format. Nevertheless, the advent of the first shaders, and with them the possibility of executing program code directly on the GPU, made it possible to use graphics accelerators for problems unrelated to graphics processing, i.e., for so-called General-Purpose computing on GPUs (GPGPU) [108]. Subsequently, GPUs became fully programmable devices. Branching operations and loops appeared, and in 2003 support for 32-bit precision computations was added. As GPUs grew more sophisticated, more advanced software development tools for them appeared as well. The chapter is organized as follows. Section 4.1 discusses the general principles and specific features of GPU programming. Section 4.2 presents a parallel program in which the host splits 37

the search tree into subtrees processed by independent GPU threads, as well as two ways of traversing the tree on the device: recursive and iterative. Section 4.3 proposes ways of optimizing the program code: using shared memory and reducing branch divergence. Section 4.4 reports the results of computational experiments and assesses the efficiency of using GPUs for the problem under consideration.

4.1

General Principles and Features of GPU Programming

One of the first programming interfaces worth mentioning is Microsoft Direct3D, which was the first to support shaders (the HLSL language). For GPU programming, the SIMD (Single Instruction, Multiple Data) approach was proposed. This approach assumes that a group of processors operating in parallel act on different data, but at any given moment all of them must execute the same instruction. An important milestone in the development of general-purpose GPU computing was the introduction in November 2006 of the Compute Unified Device Architecture (CUDA), a hardware and software platform developed by NVIDIA and released together with the GeForce 8800 GTX graphics card [109]. This technology is based on an extension of the C language and provides access to the instruction set of the graphics accelerator and control over its memory for parallel computing. Today, CUDA is one of the most popular programming technologies for computing systems equipped with GPUs [110]. The operation of a GPU can be described by the stream computing model. A stream (in the data sense) is a sequence of elements of the same type, each of which is processed independently of the others. A stream is processed by compute kernels. The body of a kernel contains the computations for a single element; it is invoked many times (possibly in parallel), once for each element. The CUDA programming model is also based on the stream computing model. An application includes code executed on the host, which implements only the sequential parts of the algorithm, prepares data and copies them to the device, sets the parameters of the kernel, and launches it. A kernel is a special function (kernel function), hundreds or thousands of instances of which run on the device in parallel, each processing its own part of the input data [111]. In effect, the GPU acts as a specialized coprocessor to the CPU, with its own memory and the ability to execute a large number of threads in parallel. Threads in CUDA are organized in a two-level hierarchy. At the upper level, a grid is a one-, two-, or three-dimensional array of blocks. In early versions of CUDA, the third dimension of block indices always had to equal 1. At the lower level, a block is a one-, two-, or three-dimensional array of threads. All threads run in parallel; inside a thread, information about the block and grid it belongs to is available through the built-in CUDA variables blockIdx (the index of the block in the grid) and threadIdx (the index of the thread in the block). Each of these variables is a three-component integer vector. The dimensions of the grid and of a block are available through the built-in variables gridDim and blockDim [110, 112]. Threads are divided into groups of 32 called warps. In fact, only threads within one warp are executed physically simultaneously, and they execute the same instructions in parallel according to the SIMD principle. Although this approach is very similar to the SIMD model, there are fundamental differences, and NVIDIA uses the term SIMT (Single Instruction, Multiple Threads). The thread scheduler periodically switches execution from one warp to another, so threads of different warps may be at different points of program execution. Warps are managed transparently by the GPU itself. For developing GPU programs, NVIDIA released its own extensions of the C language and the NVCC compiler for building such programs, and introduced the new file extension *.cu for files containing CUDA source code. The C language extensions comprise specifiers for functions and variables, new built-in data types, built-in variables (inside a kernel), and a directive for 38

launching a kernel from C code. CUDA offers two Application Programming Interfaces (APIs): the low-level driver API and the high-level runtime API. The runtime API is implemented on top of the driver API; it is less flexible but more convenient for writing programs, since it does not require a large number of initialization operations. In general, CUDA is used according to the following scheme: • allocate memory on the GPU; • copy data from CPU memory to the allocated GPU memory; • launch the kernel function (or several kernel functions in succession); • copy the results back to CPU memory; • free the allocated GPU memory. The capabilities of a CUDA GPU are characterized by its Compute Capability (CC), whose major number corresponds to the architecture version and whose minor number corresponds to incremental changes within the architecture. For example, GPUs with CC 1.0 do not support double-precision computations, GPUs with CC below 2.0 do not support recursion in kernel functions, and so on. Today, applications based on the CUDA architecture can be developed in various languages and APIs, including C, C++, Fortran, OpenCL, and DirectCompute. The CUDA architecture is actively developed and supported by NVIDIA. Its drawbacks are that it is proprietary (the technology is developed and supported by NVIDIA) and cannot be used to program graphics cards of other vendors. Alternative tools for developing GPGPU applications include the Open Computing Language (OpenCL) and Microsoft DirectX 11 DirectCompute [113, 114]. OpenCL is a framework for writing programs that perform parallel computations on various GPUs and CPUs [115]. The OpenCL framework comprises a programming language based on the C99 standard and an application programming interface (API). OpenCL provides both task-level and data-level parallelism and is an implementation of the GPGPU approach [116]. OpenCL is a fully open standard developed and maintained by the non-profit Khronos Group consortium, whose members include a number of major companies such as Apple, AMD, Intel, NVIDIA, Sun Microsystems, and Sony Computer Entertainment. Microsoft DirectCompute is an application programming interface that is part of the Microsoft DirectX API collection and is intended for computers running Microsoft Windows operating systems [117]. DirectCompute is designed for general-purpose computing on GPUs and is thus an implementation of the GPGPU concept. DirectCompute was originally released as part of DirectX 11, but later became available for DirectX 10 and DirectX 10.1 as well. In this book, NVIDIA CUDA is used to program GPGPU applications.

4.2

Development of the Parallel Program

To split the original search tree into subtrees for parallel processing, we use the following strategy, which is largely similar to the master–worker paradigm presented in Section 3.3. Here the CPU plays the role of the master, while the workers run in parallel on the GPU. The granularity parameter G is used to vary the degree of parallelism: the part of the tree above this level is processed sequentially on the CPU, while the subtrees below it are turned into parallel tasks and processed by separate threads on the GPU. The depth-first traversal of all tree nodes Ni = {ni,1 , ni,2 , . . . , ni,k }, i = G + 1, I, k = Q 1, Ki , Ki = il=1 Jl is performed by independent threads on the GPU. The total number of such Q threads is Nthreads = KG = G l=1 Jl , (K0 = 1), 0 ≤ G ≤ I − 1. In particular, Fig. 4.1 shows how the original search tree is split into subtrees for the simple CES example of Fig. 1.1. 39

Figure 4.1: Splitting the search tree into subtrees for parallel processing on the GPU For the granularity parameter G = 2, levels 0 to 2 are processed on the CPU; the initial fragments of CES variants are then transferred to the GPU, and levels 3 and 4 are processed in parallel on the GPU. The CPU does not check constraints (1.4) and (1.5), does not compute the objective function (1.1), and does not check the upper and lower bounds of the branch-and-bound method: all these operations are performed on the GPU. Host (CPU) The host (see Listing 4.1) starts by calling the procedure ReadInputData(), which reads the input data from a file (line 2). The variable inData (input data) is a pointer to a structure whose fields are variables and arrays of various dimensions holding all the necessary input data. These data are the same for all threads and do not change while the problem is being solved, i.e., they are constants. The function ThreadsNumber() computes, from the input data inData and the granularity parameter G, the number of parallel threads numThreads required to process all subtrees (line 4). This value is then used when calling the procedure PrepOperationalData(), which prepares the working data (line 6). The variable oprData (working data) is a pointer to a structure whose fields are variables and arrays of various dimensions whose values are modified independently by all threads while the problem is being solved. Therefore, one of the array indices is the thread identifier threadID. Both procedures, ReadInputData() and PrepOperationalData(), initialize the variables and perform all memory allocations needed for the arrays. Next, the procedure EnumerateHost() is called at the first level of the search tree (line 10). Using the recursive procedure EnumerateHost(), the host performs a sequential depth-first traversal of the tree down to the level given by the granularity parameter G, forming the initial fragments of CES variants and storing them in the array W (lines 33–42). W is an array of numThreads vectors of length I each. The elements of each vector hold the selected main size of the unit for each stage of the CES. Note that the host does not check constraints (1.4) and (1.5) and does not compute the upper and lower bounds of the objective function (1.1). All these operations are performed on the GPU.

40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

main () { /* prepare the input data */ ReadInputData ( inData ) ; /* d e t e r m i n e the r e q u i r e d number of threads */ numThreads = ThreadsNumber ( inData , G ) ; /* prepare the working data */ P r e p O p e r a t i o n a l D a t a ( oprData , inData , numThreads ) ; /* t r a v e r s e the tree to split it into s u b t r e e s */ /* and build the initial f r a g m e n t s of the CES */ threadID = 0; EnumerateHost (1) ; /* t r a n s f e r the input and working data to the device */ c u d a M e m c p y H o s t T o D e v i c e ( inData ) ; c u d a M e m c p y H o s t T o D e v i c e ( oprData ) ; c u d a M e m c p y H o s t T o D e v i c e (W , G ) ; /* d e t e r m i n e the kernel launch p a r a m e t e r s */ blocksPerGrid = numThreads / M A X _ T H R E A D S _ P E R _ B L O C K + 1; th re ad s Pe rB lo ck = M A X _ T H R E A D S _ P E R _ B L O C K ; /* launch the kernel f u n c t i o n on the device */ FindSolution < < < blocksPerGrid , threadsPerBlock > > >( G ) ; /* s y n c h r o n i z e with the device */ c u d a D e v i c e S y n c h r o n i z e () ; /* copy the search results from the device to the host */ cudaMemcpyDevicetoHost (W); c u d a M e m c p y D e v i c e t o H o s t ( Cost ) ; c u d a M e m c p y D e v i c e t o H o s t ( minCost ) ; /* find the global optimal s o l u t i o n */ for ( n = 0; n < numThreads ; n ++) { if ( Cost [ n ] <= minCost ) { Wopt = W [ n ]; minCost = Cost [ n ]; } } } EnumerateHost ( i ) { if ( i <= G ) { /* current search depth does not exceed G */ for ( j =1; j <= J [ i ]; j ++) { Wloc [ i ] = X [ i ][ j ]; EnumerateHost ( i + 1) ; } } /* r e c u r s i v e search */ else { /* an initial f r a g m e n t of G stages has been built ( i == G + 1) */ /* store the initial f r a g m e n t of the CES */ W [ threadID ] = Wloc ; threadID ++; } }

Listing 4.1: Pseudocode of the host (CPU) algorithm

Once the initial fragments of the CES have been formed, they are transferred to the device for further computation (lines 12–14). Then the required number of blocks in the grid, blocksPerGrid (line 16), and the number of threads per block, threadsPerBlock (line 17), are computed, and the kernel function FindSolution() is launched (line 19). When the required number of blocks is computed (line 16), the result of the division may be non-integer and is automatically truncated to an integer; therefore, 1 is added to the quotient to round it up. In CUDA, a kernel function is launched asynchronously: after the kernel has been launched on the device, control returns to the host immediately (the host does not wait for the kernel to finish). Therefore, the host must call cudaDeviceSynchronize() (line 21) to pause and wait until all GPU threads have finished (i.e., synchronize with the device). After the device has finished, the results are transferred to the host (lines 23–25). Here Cost is a one-dimensional array of numThreads elements used to store the local costs of the CES variants found by each thread. The smallest value in the array Cost of local costs is then found (lines 27–30) and stored in the variable minCost. This is the global optimal cost of the best CES variant. The optimal CES variant itself is stored in the array Wopt with I elements, each holding the main size of the unit for the corresponding stage.

41

Device (GPU) Two approaches to traversing the search tree were implemented: recursive (see Listing 4.2) and iterative (see Listing 4.3). The recursive approach is convenient for traversing tree structures and can be used on NVIDIA GPUs with Compute Capability 2.0 and higher. The iterative approach can be used on older NVIDIA devices that do not support recursion (Compute Capability below 2.0). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

_ _ g l o b a l _ _ void FindSolution ( G ) { threadID = blockDim . x * blockIdx . x + threadIdx . x ; /* if the i d e n t i f i e r is less than the number of threads */ if ( threadID < numThreads ) { En um er a te De vi c e ( threadID , G +1) ; } /* start the tree t r a v e r s a l */ } _ _ d e v i c e _ _ E nu me r at eD ev i ce ( threadID , i ) { for ( j =1; j <= J [ i ]; j ++) { /* build a c o m p l e t e CES variant */ Wloc [ i ] = X [ i ][ j ]; if ( i < I ) { /* check the o p e r a b i l i t y and the upper bound */ if ( S ( Wloc , i ) == 0 && DefineCost ( Wloc , i ) < minCost ) { En um er a te De vi c e ( threadID , i + 1) ; } } /* r e c u r s i v e search */ else { /* the last level of the tree has been reached */ /* check the o p e r a b i l i t y ( c o n s t r a i n t 1.4) and the upper bound */ if ( S ( Wloc , I ) == 0 && DefineCost ( Wloc , I ) < minCost ) { if ( T ( Wloc ) <= Tmax ) { /* check c o n s t r a i n t 1.5 */ /* store the new local optimal s o l u t i o n */ W [ threadID ] = Wloc ; Cost [ threadID ] = DefineCost ( Wloc , I ) ; atomicMinDbl ( minCost , Cost [ threadID ]) ; } } } } }

Listing 4.2: Pseudocode of the recursive device (GPU) algorithm The kernel function FindSolution() (see Listing 4.2, line 1) determines the thread identifier threadID (line 2). Its value must be less than the number of threads used, numThreads (line 4). This check is necessary because the number of blocks in the grid is rounded up when the kernel parameters are computed (see Listing 4.1, line 16). The number of threads launched on the GPU may exceed numThreads, whereas all dynamic memory allocations for the arrays in oprData were made for exactly numThreads elements. The recursive procedure EnumerateDevice() is called for all threads whose identifier is less than numThreads (line 5). Using this procedure, each thread performs a sequential depth-first traversal of the subtrees assigned to it (from level G+1 down to the lowest level I), completing the initial fragment of the CES variant received from the host to a full CES variant. To reduce the enumeration, an algorithm based on branch-and-bound is used here as well. The set of feasible solutions is successively partitioned into subsets [68]. A lower bound on the objective function is then determined for each subset. If the lower bound is not less than the best solution found so far, the subset is discarded, and the process continues until a single solution is obtained (assuming that the optimality criterion is the minimum of the objective function). If the current tree level i is not the lowest level I (line 11), the compatibility constraint (1.4) for the joint operation of the CES equipment is checked at level i (procedure S(), line 13) for the initial fragment of the CES consisting of the levels already passed above and the current one (the CES fragment consisting of stages 1 to i). If part of the flowsheet is already inoperable, there is no point in enumerating units at the lower levels. In addition, the cost of the initial fragment of the CES variant (procedure DefineCost(), line 13) is compared with the cost minCost of the current optimal CES variant. As a result, the number of enumerated variants is reduced. The variable minCost can be read and written by all threads running on the device. Its initial value is the sum of the maximum unit costs over all stages. The algorithm descends to the next recursion level (a call of EnumerateDevice() at level i + 1, line 14) only if the initial fragment of the CES variant consisting of stages 1 to i is, first, 42

operable and, second, cheaper than the best CES variant found so far; otherwise, it backtracks to the previous level of the tree. Thus, when the last level of the tree is reached (i == I, line 15), a complete operable CES variant (one satisfying constraint (1.4)) has been formed. Its cost (the result of DefineCost()) is then compared with the cost minCost of the current best CES variant (line 17). If a better solution has been found, the processing time constraint (1.5) is checked (procedure T(), line 18). If constraint (1.5) is satisfied, a new optimal CES variant has been found; it replaces the current optimal solution W (line 20), and its cost is stored in minCost using the atomic function atomicMinDbl() (line 21). Atomicity guarantees that the operation is performed correctly when many threads attempt to execute it simultaneously. A well-known way of traversing a tree without explicit recursion is the iterative use of a stack [118]. A stack is a data structure representing a list of elements organized according to the LIFO (Last-In, First-Out) principle. The pseudocode of the algorithm implementing this approach is shown in Listing 4.3. The kernel function FindSolution() is identical to the one in Listing 4.2 (lines 1–5). 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

... _ _ d e v i c e _ _ E nu me r at eD ev i ce ( threadID , i ) { /* declare a s t a n d a r d Last - In , First - Out ( LIFO ) stack */ stack [ MA X_STACK_ SIZE ]; stackTop = 0; /* define the s t a n d a r d stack o p e r a t i o n s */ # define push ( x ) stack [ stackTop ++] = ( x ) # define pop () stack [ - - stackTop ] # define empty ( stackTop == 0) /* push the first node of level i onto the stack */ push (1) ; /* repeat while the stack is not empty */ while (! empty ) { j = pop () ; /* pop the next node from the stack */ /* if there are no more nodes at this level , go one level up */ if ( j == J [ i ] + 1) {i - -; } else { /* o t h e r w i s e push the next node of the current level */ push ( j + 1) ; /* build a c o m p l e t e CES variant */ Wloc [ i ] = X [ i ][ j ]; /* check the o p e r a b i l i t y and the cost of the CES variant */ if ( S ( Wloc , i ) == 0 && DefineCost ( Wloc , i ) < minCost ) { if ( i < I ) { /* the last level has not been reached */ i ++; /* go to the next level of the tree */ /* push the first node of the next level */ push (1) ; } else { /* the last level has been reached ( i == I ) */ /* check the o p t i m a l i t y c r i t e r i o n */ if ( DefineCost ( Wloc , I ) < minCost ) { /* check the p r o c e s s i n g time c o n s t r a i n t (1.5) */ if ( T ( Wloc ) <= Tmax ) { /* store the new local optimal s o l u t i o n */ W [ threadID ] = Wloc ; Cost [ threadID ] = DefineCost ( Wloc , I ) ; atomicMinDbl ( minCost , Cost [ threadID ]) ; } } } } } } }

Listing 4.3: Pseudocode of the iterative device (GPU) algorithm The stack is implemented as a one-dimensional static array stack of length MAX_STACK_SIZE (line 10). The stack size (the value of MAX_STACK_SIZE) must be at least the number of tree levels: MAX_STACK_SIZE >= I. Data can be added to or removed from the stack only at its top. The variable stackTop points to the top of the stack. Next, the standard elementary stack operations are defined. The push operation adds a new element (line 13), the pop operation removes the element at the top of the stack (line 14), and the empty operation checks whether the stack is empty (line 15). 43

Initially, the stack is empty, and the first tree node at level G + 1 (the first unit variant for stage G + 1 of the CES) is pushed onto it (line 17). Then the tree nodes are enumerated iteratively until the stack becomes empty (line 19). At each iteration of the loop, the most recently pushed node is popped from the stack (line 20). If there are no more nodes at the current level, the algorithm moves one level up (line 22). Otherwise, the next node of the current level is pushed onto the stack (line 25). If the current initial fragment of the CES consisting of stages 1 to i satisfies the compatibility condition (1.4) and its cost is lower than the cost of the best CES variant found so far (line 29), the unit variant at the current level is appended to the preceding operable initial fragment of the CES consisting of stages 1 to i-1 (line 27). Otherwise, the algorithm does not descend to deeper levels of the tree and proceeds to the next loop iteration. If the current level is not the last one (line 30), the algorithm moves to the next level (line 31) and pushes its first node onto the stack (line 33). If the last level of the tree has been reached (line 34), the algorithm proceeds in the same way as the recursive one (see Listing 4.2). The cost of the operable CES variant (the result of DefineCost()) is compared with the cost minCost of the current best CES variant (line 36). If a better solution has been found, the processing time constraint (1.5) is checked (procedure T(), line 38). If constraint (1.5) is satisfied, a new optimal CES variant has been found; it replaces the current optimal solution (line 40), and its cost is stored (lines 41–42). For the computational experiments, the algorithm variants presented above were implemented in the C programming language. They use the mathematical model of CES operation (including, among other things, the expressions for checking the constraints and computing the optimality criterion) described in detail in [58]. The amount of input data required for the CES equipment selection problem is small (tens of kilobytes), but the data have a rather complex structure. In particular, in [54, 55, 119] the model of CES operation was implemented using an objectoriented approach (the C++ programming language) and the Standard Template Library (STL), in particular the std::vector template whose elements were other vectors. For example, a CES is a vector of process routes, a process route is a vector of processing stages, a processing stage is a vector of operating cycles, an operating cycle is a vector of process operations, and so on. However, at the time of this study, CUDA had no full-fledged implementation of STL containers that could be used in code executed on the GPU. For this reason, the program code was completely rewritten without STL. Instead of numerous nested std::vector objects, arrays (about 50 of them) of various dimensions (from 1 to 5) were used. To make the arrays occupy contiguous regions of memory (flat memory), multidimensional arrays were converted into one-dimensional ones (see Listing 4.4). This guarantees that the data occupy a single contiguous memory region, which is convenient when transferring data with memcpy-like functions. Moreover, since the data are not fragmented in memory, the number of cache misses decreases, which considerably increases the program’s performance. On the test platform described below, the new version of the sequential CPU program (without nested vectors) turned out to be 8–9 times faster than the previous STL-based version. In what follows, all computational experiments (performance comparisons, speedup calculations, etc.) are carried out with the new version of the program.

4.3

Optimization of the Program Code

To increase the performance of the program, two optimization approaches were applied: efficient use of shared memory and reduction of branching in the kernel function. Physically, GPU memory can be divided into Dynamic Random-Access Memory (DRAM) and memory located directly on the chips of the streaming multiprocessors (SMs). However, the classification of memory is not limited to its physical location.

44

... typeName * A ; /* A [ i ][ j ] , a two - d i m e n s i o n a l array of size I * J */ typeName * B ; /* B [ i ][ j ][ k ] , a three - d i m e n s i o n a l array of size I * J * K */ /* a l l o c a t e memory */ A = malloc ( I * J * sizeof ( typeName ) ) ; B = malloc ( I * J * K * sizeof ( typeName ) ) ; /* i n i t i a l i z e the arrays */ for ( i = 0; i < I ; i ++) { for ( j = 0; j < J ; j ++) { A[i * J + j] = i + j; for ( k = 0; k < K ; k ++) { B[i * J * K + j * K + k] = i * j * k ;}}} ...

Listing 4.4: Accessing multidimensional arrays as one-dimensional ones

The memory hierarchy of NVIDIA GPUs of the Kepler and Fermi microarchitectures is shown in Fig. 4.2 [120, 121].

Figure 4.2: Memory hierarchy of NVIDIA GPUs of the Kepler and Fermi microarchitectures (based on [120, 121]) Global memory is physically implemented as DRAM chips soldered directly onto the graphics board. It is the largest memory on the card, measured in gigabytes. However, since it is located outside the multiprocessors, it is the slowest type of memory [110]. In the Fermi and Kepler microarchitectures, global memory is cached by the L2 and L1 caches, but it nevertheless remains the slowest (with a latency of 400 to 600 clock cycles without caching). Global memory can be read and written by all threads of the grid. Optimizing global memory accesses usually has a very large effect on program performance. In addition, each streaming multiprocessor has a small amount (64 KB) of on-chip memory. In the Fermi microarchitecture, the programmer can configure this memory either as 48 KB of shared memory and 16 KB of L1 cache, or as 16 KB of shared memory and 48 KB of L1 cache. The Kepler microarchitecture also allows this memory to be split in half: 32 KB of shared memory and 32 KB of L1 cache. The lifetime of data in shared memory coincides with the lifetime of the block. In addition to the L1 cache, the Kepler microarchitecture has a 48 KB 45

read-only data cache. The L2 cache, shared by all SMs, is the main point at which data and instruction read/write requests between the streaming multiprocessors and DRAM are unified. Constant memory physically resides in DRAM (like global memory), but is cached for reading on all devices. Constant memory is read-only for the GPU and readable and writable for the CPU. In the first implementation of the program, all input data inData were loaded into constant memory and the working data oprData into global memory. Shared memory is the fastest, since it is located directly on the SM chip. Therefore, in the first optimization (optimization Mo1 ), the kernel function FindSolution() was modified (see Listing 4.5) so as to move the local array Wloc (line 7) and, for the iterative algorithm, also the array stack (line 9) from global memory to shared memory using the procedure SendDataByMemoryPtr() [122]. 3 4 5 6 7 8 9 10 11 12 13 14 15

... threadID = blockDim . x * blockIdx . x + threadIdx . x ; _ _ s h a r e d _ _ sharedMemory []; /* copy data to shared memory */ S e n d D a t a B y M e m o r y P t r ( Wloc , sharedMemory ) ; # ifdef ITERATIVE S e n d D a t a B y M e m o r y P t r ( stack , sharedMemory ) ; # endif /* if t h r e a d I D is less than the number of threads */ if ( threadID < numThreads ) { /* start the tree t r a v e r s a l */ En um er a te De vi c e ( threadID , G +1) ; } ...

Listing 4.5: Pseudocode of the optimized algorithm: optimization Mo1 Both arrays, Wloc and stack, are used intensively during the traversal of the decision tree. In the second variant (optimization Mo2 ), the input data inData, whose size is several tens of kilobytes, were moved from constant memory to shared memory (see Listing 4.6, line 7). For this purpose, the procedure SendDataByMemoryPtr() is used; it runs on the GPU and transfers data to the given memory address sharedMemory. 3 4 5 6 7 8 9 10 11 12

... threadID = blockDim . x * blockIdx . x + threadIdx . x ; _ _ s h a r e d _ _ sharedMemory []; /* copy data to shared memory */ S e n d D a t a B y M e m o r y P t r ( inData , sharedMemory ) ; /* if t h r e a d I D is less than the number of threads */ if ( threadID < numThreads ) { /* start the tree t r a v e r s a l */ En um er a te De vi c e ( threadID , G +1) ; } ...

Listing 4.6: Pseudocode of the optimized algorithm: optimization Mo2 The amount of extremely fast shared memory sharedMemory[] is specified by the programmer as a kernel launch parameter, but cannot exceed 48 KB because of GPU hardware limitations. The threads, grouped into blocks by the programmer, are executed on the GPU. For each SM separately, the warp scheduler divides the blocks into warps, and all threads of a warp always belong to the same block. It is within a warp that threads are executed physically simultaneously. In the CUDA implementation used, a warp consists of 32 threads. An SM executes warps one after another on its streaming processors, and switching between warps incurs minimal overhead. All threads execute the same instruction but process different data. Depending on the program logic and on the data of each thread, different threads may be in different states at the same moment. If the control flow of different GPU threads follows different execution paths, so-called branch divergence occurs. A warp executes all branches of the kernel code serially, disabling the threads whose code 46

is not in the current branch [123, 124]. The fewer strongly divergent branches (i.e., branches executed differently depending on the arguments), the higher the performance. The most common cause of divergence is the use of an if-then-else statement in the kernel code. If the condition evaluates to true for some threads of a warp and to false for the others, the execution path of the control flow branches. The problem is that the then and else branches are executed serially rather than in parallel. First, all instructions of the then branch are executed, while all threads for which the condition is false are suspended. Then the instructions of the else branch are executed, and all threads for which the condition is true are suspended in turn. This causes a considerable loss of performance, which cannot be compensated for by a higher degree of parallelism [123, 125]. Besides conditional statements, this problem can be caused by irregular loops, i.e., loops whose number of iterations is not known in advance, such as loops with a precondition or a postcondition (while, do-while) [126]. The number of iterations of such loops may differ from thread to thread, so that some threads stay idle waiting for the loop to finish in other threads [127]. The conditions that cause the control flow to branch often depend on the data and cannot be determined at compile time. In this book, we attempt to reduce the effect of branch divergence by removing from the kernel function the expressions that check the operability of the CES variant and the upper bound, i.e., the comparison of the cost of the current CES variant with the cost of the best variant found so far (see Listing 4.2, line 13, and Listing 4.3, line 29). All conditions are now checked only at the lowest level of the tree (i==I), so the execution paths of the control flow branch less. Removing the checks at each level of the tree traversal effectively turns the branch-and-bound scheme into exhaustive search [9]. The names of the algorithm versions optimized for divergence reduction carry the suffix “dr” (Mo1dr and Mo2dr).

4.4

Results of Computational Experiments

For the computational experiments, a program was developed in C using CUDA [109]. The experiments were carried out on a computing system with an Intel Xeon processor (E5-1620 v2, 4 cores with Hyper-Threading, 10 MB Smart Cache, 3.7 GHz) and 16 GB of RAM, equipped with an NVIDIA Tesla K20c GPU. GPU specifications: 13 streaming multiprocessors (SMs) with 192 CUDA cores each (2496 CUDA cores in total), 5 GB of global memory, and 48 KB of shared memory per SM. The software comprised Ubuntu 14.04.1, NVIDIA driver 340.29, CUDA 6.5, and the GNU C++ compiler 4.8.2. The test case was the CES equipment selection problem described in Section 1.3 (see Fig. 1.4), with 3 standard equipment sizes available for each of the 16 processing stages (316 = 43 046 721 CES variants in total). The program was run with values of the granularity parameter G from 0 to 10. The value G = 0 corresponds to the case where the host forms no initial fragments of CES variants and the whole search tree is traversed by a single GPU thread (Nthreads = K0 = 1), i.e., to sequential execution of the algorithm on the GPU. This value serves as the baseline for assessing the effect of parallelization. For G > 10, too many threads are created on the GPU, and the working data oprData require more memory than the available DRAM (an out-of-memory error occurs). In optimization Mo1, shared memory is used to store the array Wloc in the recursive approach, and to store Wloc and stack in the iterative approach. For granularity values greater than 6 (recursive approach) and greater than 5 (iterative approach), a shared memory shortage error occurs. The results of the Mo1 and Mo2 optimizations are shown in Fig. 4.3 (the vertical axis is logarithmic). The total run time of a program on a heterogeneous computing system consists of the run 47

100000

Recursive Recursive Mo1 Recursive Mo2 Iterative Iterative Mo1 Iterative Mo2

Run time, s

10000 1000 100 10 1 0.1

0

1

2

3

4

5

6

Granularity parameter G

7

8

9

10

Figure 4.3: Shared memory optimization: run time of the parallel program versus granularity time of the sequential part on the CPU and the run time of the parallel part on the GPU. The experiments showed that the sequential part on the CPU accounts for a negligible fraction of the total run time (no more than 0.1%). For G > 10, global memory runs out for all versions of the program (for G > 10, at least 311 = 177 147 threads are required). For G > 6 for the “Recursive Mo1 ” version and for G > 5 for the “Iterative Mo1 ” version, shared memory runs out. Optimization Mo1 had no significant effect on the total run time (the average change in run time was less than 1%). Optimization Mo2 reduced the run time by 48% on average. It should also be noted that the recursive approach turned out to be faster than the iterative one by 120% on average, i.e., other things being equal, the recursive approach is on average 2.2 times faster than the iterative one, which uses a stack. The results of reducing the amount of branching are shown in Fig. 4.4 (the vertical axis is logarithmic). In this case, the branching function was removed as described above, and a series of computational experiments similar to that for the memory optimization was carried out. In general, the pattern of the results remained the same as before: optimization Mo2dr turned out to be the most effective. At the same time, removing the branching function from the algorithm had markedly different effects on the recursive and iterative approaches. For the recursive algorithm, the run time increased slightly (by no more than 5–6%) compared with the original version. For the iterative approach, the run time decreased by 55–58% on average. The experimental data suggest that, in terms of run time reduction, optimization Mo2 is the most effective for the recursive approach, and Mo2dr for the iterative approach. For the test case used, the fastest recursive version, Mo2, is on average 5% faster than the fastest iterative version, Mo2dr. A number of tests for larger problems were also carried out to measure the speedup relative to a fully sequential CPU program. The experiments used the same test CES, with 3 (316 = 43 046 721 CES variants in total), 4 (416 = 4 294 967 296 variants), and 5 (516 = 152 587 890 625 CES variants) standard equipment sizes available for each stage. The results for the recursive and iterative approaches are shown in Fig. 4.5. 48

100000

Recursive dr Recursive Mo1dr Recursive Mo2dr Iterative dr Iterative Mo1dr Iterative Mo2dr

10000

Run time, s

1000 100 10 1 0.1

0

1

2

3

4

5

6

Granularity parameter G

7

8

9

10

Figure 4.4: Branching reduction: run time of the parallel program versus granularity For the recursive implementation, the speedup ranges from 2.66 to 5.79, and for the iterative implementation, from 1.67 to 4.75. In all computational experiments, the recursive approach was faster than the iterative one. Although the program as a whole accelerates the solution of the problem considered in this book, a direct implementation of branch-and-bound on the GPU does not yield a large speedup (by hundreds or thousands of times). The reason is that the traversal of the search tree involves a great deal of branching in the program code, which considerably reduces the efficiency of the GPU. Similar problems are described, for example, in [128], where the authors report the results of implementing minimax search using CUDA.

49

6

Recursive Mo2 Iterative Mo2dr

5

Speedup

4 3

2

1 0

316

416

516

Number of CES variants (problem size)

Figure 4.5: Speedup for problems of different sizes

50

Chapter 5

Parallel Hybrid Metaheuristic Algorithm on GPUs As shown in Chapters 2–4, branch-and-bound is guaranteed to find the globally optimal solution, but its run time grows rapidly with the problem size. Despite the optimization of the sequential algorithm and its parallelization for shared-memory and distributed-memory systems, solving large problems still takes tens of hours, and a direct implementation of branch-and-bound on the GPU does not yield a substantial speedup (Chapter 4). In this chapter, we therefore propose to solve the problem with metaheuristic algorithms, which do not guarantee finding the global optimum but produce near-optimal solutions in acceptable time. The algorithm proposed in this book combines two metaheuristics: Ant Colony Optimization (ACO) and Simulated Annealing (SA) [129]. In the proposed approach, ACO is used to find a feasible initial solution for SA [130]. This division of roles is due to the specific features of the problem. For classical combinatorial optimization problems such as the travelling salesman problem, a randomly generated variant can serve as the initial solution. For the CES equipment selection problem, this approach is not applicable: the initial solution must satisfy both the compatibility constraint (1.4) and the processing time constraint (1.5), and the probability that a randomly chosen variant is feasible is low [131]. Therefore, ACO is used for a targeted search for a feasible CES variant, and SA is used to improve it subsequently with respect to the capital cost of the equipment. Both algorithms are stochastic and well suited to the massively parallel architecture of the GPU: each thread can search for a solution independently, and the results of all threads are then compared. This makes it possible to explore a much larger part of the search space in the same time than with sequential execution. The chapter is organized as follows. Section 5.1 gives a brief overview of metaheuristic algorithms for combinatorial optimization problems. Sections 5.2 and 5.3 describe Simulated Annealing and Ant Colony Optimization, together with the ways of constructing a candidate solution and defining the heuristic information and the fitness function proposed for the problem under consideration. Section 5.4 presents the implementation of the hybrid algorithm (ACO+SA) on the GPU, and Section 5.5 reports the results of computational experiments and a comparison with branch-and-bound.

5.1

A Brief Overview of Metaheuristic Algorithms for Combinatorial Optimization

Metaheuristic algorithms combine basic heuristic methods within higher-level algorithmic frameworks aimed at an efficient exploration of the search space. The main difficulty lies in adapting general metaheuristic frameworks to specific problems. The fundamental properties of meta-

51

heuristics include the following [132]: • metaheuristics are strategies that guide the search process; • the goal of a metaheuristic is to explore the search space efficiently in order to find near-optimal solutions; • metaheuristic algorithms are approximate and, as a rule, nondeterministic; • metaheuristic algorithms may include mechanisms for avoiding being trapped in a confined region of the search space; • metaheuristics are not necessarily problem-specific and can be described at an abstract level; • metaheuristics use domain-specific knowledge in the form of heuristics controlled by a higher-level strategy; • modern metaheuristics use search experience stored in memory to guide the search. Each metaheuristic algorithm has its own domain-dependent features; nevertheless, all of them share a number of invariable basic components, which can be represented at an abstract level as a generic scheme. This generic scheme is as follows [132, 133]. 1. Initialization: finding an initial solution. 2. Forming one or more local neighborhoods of the initial solution. 3. Selecting the most promising neighborhood. 4. Selecting a candidate in the chosen neighborhood using an evaluation function (fitness function) that depends on such parameters of two solutions as the objective function value, penalties for violating certain constraints, etc. The result of the selection is the best solution with respect to this criterion. 5. Stopping according to one of the criteria, which may be the computation time, the number of iterations, or the rate of improvement of the solutions. Any classification of metaheuristic methods is somewhat arbitrary, since typical groups of methods are based on similar ideas and there is a persistent trend towards hybrid algorithms [23, 134]. Four groups of metaheuristic algorithms are distinguished in [135]: evolution-based algorithms, swarm intelligence-based algorithms, physics-based algorithms, and multi-start algorithms. Evolution-based algorithms Genetic algorithms (with binary and real-valued encoding); methods that imitate the immune systems of organisms (artificial immune systems); scatter search; the covariance matrix adaptation evolution strategy; variable mesh optimization; differential evolution; weed colonization; and cuckoo search. Swarm intelligence-based algorithms Particle swarm optimization; ant colony optimization; bacterial foraging optimization; bee colony algorithms, including the artificial bee colony algorithm and the bees algorithm; fish school search; the bat-inspired algorithm; the firefly algorithm; glowworm swarm optimization; and the shuffled frog-leaping algorithm. 52

Physics-based algorithms Central force optimization; simulated annealing; adaptive simulated annealing; deterministic annealing; harmony search; and the electromagnetism-like mechanism. Multi-start algorithms The greedy randomized adaptive search procedure (GRASP) and tabu search.

5.2

Simulated Annealing (SA)

Simulated Annealing (SA) [136], an optimization technique that uses a controlled random search to avoid getting trapped in local optima, was chosen as the algorithm for minimizing the objective function (the equipment cost). The idea of the method arose from observing the processes that take place in a metal as it cools after being heated to its melting point [137, 138]. Simulated annealing is a stochastic search method in which, at each step, the current solution is replaced by another solution randomly chosen from its neighborhood. A solution that improves the objective function is always accepted, while a worsening solution is accepted with a probability that depends on a control parameter called the temperature. The temperature is gradually lowered according to a cooling schedule, so that worsening solutions are often accepted at the beginning of the search and are practically never accepted towards its end. Simulated annealing can be used to search for solutions of complex linear and nonlinear problems and is able to escape from local minima. Its advantages include the ability to search for the global optimum, its generality, and its relative ease of implementation. The search strategy allows the algorithm to jump over local minima of the objective function and reach the basin of attraction of the global minimum [18, 139]. This basin is reached only with a certain probability; therefore, simulated annealing does not guarantee that the global minimum will be found. However, with a proper choice of the “temperature” trajectory (the parameter of the method), the algorithm not only improves the initial approximation but also comes substantially closer to the global optimum. The general scheme of simulated annealing is shown in Listing 5.1. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

S i m u l a t e d A n n e a l i n g () { /* i n i t i a l i z e the t e m p e r a t u r e and the initial s o l u t i o n */ t = Tinit ; W = Winit ; while ( t > Tfinal ) { /* until the final t e m p e r a t u r e is reached */ for ( l = 0; l < Lmax ; l ++) { /* repeat Lmax times */ Wcand = C o n s t r u c t S o l u t i o n ( W ) ; /* c o n s t r u c t a new c a n d i d a t e s o l u t i o n */ deltaCost = DefineCost ( Wcand ) - DefineCost ( W ) ; if ( deltaCost < 0) { /* if the c a n d i d a t e is better */ W = Wcand ; } /* accept the c a n d i d a t e s o l u t i o n */ else { /* o t h e r w i s e */ p = exp ( - deltaCost / t ) ; /* compute the p r o b a b i l i t y */ r = rand (0 , 1) ; /* g e n e r a t e a random number */ if ( p > r ) { /* if the p r o b a b i l i t y exceeds the random number */ W = Wcand ; } } } /* accept the c a n d i d a t e s o l u t i o n */ t = sigma * t ; } /* lower the t e m p e r a t u r e */ }

Listing 5.1: Simulated Annealing (SA) At the start of the algorithm, the temperature t is initialized with the initial temperature Tinit, and the current solution W with a feasible initial solution Winit (line 3). The feasible initial solution is the result of the ant colony algorithm. A solution W is a vector whose length equals the number of CES stages. Each element of the vector holds the index of one of the units available for installation at the corresponding stage.

53

The cooling schedule is the law according to which the temperature changes during cooling. It is defined by the initial temperature Tinit (line 3), the final temperature Tfinal (line 4), and the cooling rule, which changes the temperature at each cooling iteration (line 15). In practice, a good cooling schedule is designed by trial and error, by trying different parameter values and observing the results. Nevertheless, certain recommendations are given in the literature [18]. The initial temperature Tinit is defined as the largest possible difference in cost between two feasible solutions [18]. For the problem under consideration, the maximum difference between solutions was computed as follows: first, the most expensive unit is installed at each CES stage, and the cost of the resulting solution is calculated. This is the cost of the most expensive solution. The initial temperature is thus the difference between the potentially most expensive and the cheapest solution. The final temperature Tfinal is defined as the smallest possible difference in cost between two neighboring solutions [18]. For the problem under consideration, the final temperature was taken to be the minimum unit cost over all CES stages. This value is of the same order as the minimum nonzero difference between the costs of two units of the same stage. The algorithm contains two nested loops. The outer loop (line 4) gradually lowers the temperature t, which reduces the probability of accepting solutions that do not improve the objective function (1.1). At each particular temperature, Lmax iterations of the algorithm must be performed; this is done by the inner loop (line 5). The number of iterations Lmax affects the speed and convergence of the algorithm: if it is too large, the run time increases; if it is too small, the solutions obtained are far from the global minimum. In general, simulated annealing generates a finite inhomogeneous Markov chain. A Markov chain is a sequence of trials in which the probability of the outcome of each trial depends only on the outcome of the preceding trial. In simulated annealing, a trial consists in choosing the next solution at each step. In practice, an infinite chain cannot be executed, so its length must be limited [18]. There is no single generally accepted approach to determining Lmax ; different researchers recommend different choices. For example, Skiena [140] recommends taking the number of iterations before each temperature reduction between 100 and 1000. As described below, a new solution in the neighborhood of the current feasible solution is found by randomly changing one randomly chosen stage of this solution. For each CES stage, a set of available equipment is given. We propose the following approach to determining Lmax : the number of repetitions at the current temperature equals the total number of units available for the whole CES. For example, consider a CES of 5 stages in which 3 standard unit sizes are available for each stage. The total number of units is 3 + 3 + 3 + 3 + 3 = 15, i.e., Lmax = 15. This gives a theoretical chance of trying every available unit at every stage of the current solution. This is not an enumeration of combinations, which has exponential complexity, because at each iteration we randomly change only one randomly chosen stage of the solution, keeping the others fixed. In P general, Lmax = Ii=1 Ji . An advantage of this approach is that Lmax adapts to the problem size: the more units per stage and the more stages, the more iterations are performed before the temperature is lowered. At each iteration of the inner loop (line 5), a candidate solution Wcand is constructed in the neighborhood of the current feasible solution (line 6). For this purpose, the procedure ConstructSolution() is used. Then the change in the objective function, deltaCost, is determined (line 7). If a better solution has been obtained (line 8), it is accepted as the current solution (line 9). Otherwise, the probability p of accepting the worse candidate solution as the current one is computed (line 11). Historically, the first annealing scheme was the so-called Boltzmann annealing [136]. The probability is computed as p = exp(−∆E/(kB · t)), where kB is the Boltzmann constant and

54

∆E is the change in the energy level. The probability formula has a certain thermodynamic meaning [140]. The probability of a transition from a high-energy state to a low-energy state is very high. Nevertheless, there is a nonzero probability of a transition to a higher-energy state, and small transitions are more likely than large ones. Moreover, the higher the temperature, the higher the probability of energy transitions. For optimization problems, this means that at the beginning of the search, when the temperature is high, non-optimal solutions are accepted as current more often, and as the system cools, the probability of accepting non-optimal solutions decreases with the temperature. This is a valuable property of simulated annealing, since it allows the search to escape from local optima and, at least theoretically, to find the global optimum. In practical problems, the change in the energy level is usually taken to be the change in the objective function, and kB = 1 [15]. Thus, the transition probability formula takes the form p = exp (-deltaCost / t) used in the algorithm (line 11). After the probability p has been computed, a random number r between 0 and 1 is generated (line 12). If p exceeds r, the candidate solution is accepted as the current solution (line 14). The temperature t is then lowered according to the chosen cooling rule (line 15). There are various ways of lowering the temperature, for example, linear, t = t − ∆t, or exponential, t = σ · t. Cooling is faster with the exponential rule, so we use it, following the recommendations of [138]. The constant is usually chosen in the range 0.8 ≤ σ ≤ 0.99 [18]; for the problem under consideration, we take σ = 0.9. Note that, at the initial step, the current solution must be initialized with a feasible initial solution (a feasible CES variant W at the temperature Tinit). For typical problems such as the travelling salesman problem, the initial solution can be chosen at random. For the problem under consideration, a random solution amounts to choosing, at random, a unit for each CES stage from the set of units available for that stage. A CES variant can be generated in this way, but it may be either inoperable (violating the compatibility constraint (1.4)) or unable to produce the required products within the given time (violating the processing time constraint (1.5)). The time needed to obtain an operable CES variant in this way is completely unpredictable, so this approach is not applicable in practice. Instead, a Constraint Satisfaction Problem (CSP) is solved, in which the solution cost (1.1) is ignored and only constraints (1.4) and (1.5) must be satisfied [131, 141]. The constraint satisfaction problem can be formulated as follows: find an operable CES variant that satisfies the compatibility constraint (1.4) and the processing time constraint (1.5). Problems of this type can be solved by Ant Colony Optimization (ACO), presented in Section 5.3. The result is a CES variant that is operable and able to produce the required products within the given time, although its cost may be far from optimal. This CES variant is used as the starting point for simulated annealing, which in turn searches for a cost-optimal CES variant.

5.3

Ant Colony Optimization (ACO)

The main idea of this approach is to model the behavior of ants searching for a path from their colony to a food source [142]. Ants are able to find the shortest path from the nest to a food source quickly and to adapt to changing conditions [143, 144]. An ant colony can be regarded as a multi-agent system in which the agents (ants) interact to achieve a common goal. The ants interact indirectly through a chemical substance, the pheromone, which they deposit along their path (trail). The pheromone concentration on a trail determines how attractive it is (when choosing a path, an ant is more likely to choose the one with a higher pheromone concentration). The pheromone evaporates over time, which allows the ants to adapt their behavior to changes in the environment. The pseudocode of the ant colony algorithm is shown in Listing 5.2. The flag isFound, which controls the search, is initialized to false (line 2). The search 55

1 2 3 4 5 6 7 8 9 10 11 12

A n t C o l o n y O p t i m i z a t i o n () { isFound = false ; iterCounter = 0; Initialize () ; /* i n i t i a l i z e the p h e r o m o n e c o n c e n t r a t i o n */ /* until a s o l u t i o n has been found */ while (! isFound && iterCounter < m a x I t e r a t i o n N u m b e r ) { foreach ( ant in swarm ) { /* for each ant in the colony */ C o n s t r u c t S o l u t i o n ( alpha , beta ) ; } /* c o n s t r u c t a new s o l u t i o n */ if ( isFound ) return ; /* stop if a s o l u t i o n has been found */ Ph er om o ne Up da t e () ; /* update the p h e r o m o n e c o n c e n t r a t i o n */ E v a p o r a t e P h e r o m o n e ( rho ) ; /* p h e r o m o n e e v a p o r a t i o n */ iterCounter ++; /* i n c r e m e n t the i t e r a t i o n counter */ }

Listing 5.2: Ant Colony Optimization (ACO)

continues until a feasible solution is found or the number of attempts exceeds the maximum number maxIterationNumber (line 5). The local iteration counter iterCounter, initialized to 0 (line 2), protects the algorithm from looping forever if no ant of the colony manages to find a feasible solution (which is quite possible for a stochastic algorithm). In this case, iterCounter exceeds maxIterationNumber, and the algorithm terminates without a result; it is then restarted with new random initial values of the pheromone concentration. Before the algorithm starts, the pheromone concentration is initialized (line 3). While constructing a solution of the problem (line 7), each ant of the colony (line 5) can move in the search tree (see Fig. 1.2) only from top to bottom. Suppose that an ant is at node r = ni,j at level i of the tree. It can then move only to a child node s = ni+1,j ′ , s ∈ Cr , at level i + 1, where Cr is the set of child nodes of node r. The complete path of an ant from level 1 to level I corresponds to a potential solution of the problem. The probability prs that an ant moves from node r to node s is determined from two types of information, the pheromone level τrs and the heuristic information ηrs , according to the random β α · η β /P α proportional rule prs = τrs rs k∈Cr (τrk · ηrk ), where Cr is the set of child nodes of node r, k is the index of a child node, and α and β are empirical parameters of the algorithm that determine the relative influence of the pheromone concentration and of the heuristic information, respectively, on the probability of choosing a node [130, 138, 142]. As a rule, the values of α and β are chosen from the interval (0, 2) [145]. For α = 0, only the heuristic information is taken into account. For β = 0, only the pheromone concentration is taken into account, so that the search degenerates: all ants follow a single suboptimal route (the path with the highest pheromone concentration, which they thereby increase further). The first version of the algorithm used the values α = 0.4 and β = 0.6, obtained empirically: the sequential algorithm was run just a few times with different values of α and β, and the values yielding the highest speed were chosen. Looking ahead, we note that these values were not the best ones. The problem of finding the best values of the algorithm parameters is considered in detail in Chapter 6. The node to move to is selected by the “roulette wheel” principle: each child node has its own sector of the wheel with an area proportional to its probability. To select a node, a random number is generated, and the sector into which it falls is determined; the node corresponding to this sector is chosen. For the problem of finding a feasible CES variant, we propose the following way of defining the heuristic information, based on the following considerations. First, it is known that a CES whose stages are equipped with larger units is more expensive, but has a larger batch size and, as a rule, produces the required amount of product faster, since fewer batches are needed. For the same required production amount, a small batch size means producing many batches over a long time, whereas a large batch size means producing few batches over a short time. For example, suppose that 1000 kg of product must be produced. If the batch size (the amount of product produced in one CES cycle) is 10 kg, then 1000/10 = 100 batches must be produced (100 CES cycles), whereas with a batch size of 100 kg only 1000/100 = 10 batches are needed. The fewer 56

the cycles (i.e., the larger the units and hence the batch size), the more likely the processing time constraint (1.5) is to be satisfied. Thus, when a unit is chosen for the next stage, preference should be given to larger units. Second, the compatibility condition (1.4) must be satisfied for all initial fragments of the CES consisting of stages 1–2, 1–2–3, . . . , 1–2–3–. . .–I. Suppose that, while constructing a solution, an ant is at stage i of the CES. Then the operability of the initial fragment of the CES formed at the preceding tree levels is checked. If the corresponding CES fragment is inoperable, its weight is zero; otherwise, it is one. Thus, the heuristic ηrs for the move from node r to node s is the sum of the main size of the unit on the tree edge (r, s) and a value (0 or 1) characterizing the operability of the initial fragment of the CES. The second term is the same for all child nodes of node r. After an ant of the colony has reached the last level I of the search tree, the pheromone level on all edges of the ant’s path must be updated. The amount of pheromone deposited by an ant depends on the quality of the solution obtained (how close the solution is to the goal). Each ant deposits pheromone on the edges of its own path (line 9): for every edge (r, s) traversed by the mth ant, τrs = τrs + Q/Lm , m = 1, M , where Q is a constant, Lm is the fitness value of the solution constructed by the mth ant, and M is the number of ants in the colony. This rule has the form of the pheromone deposit rule of classical ant colony algorithms for the travelling salesman problem, where Lm is the length of the tour of the mth ant and Q is chosen to be of the same order as the length of the optimal tour [142, 143]. Here Lm is the value of the fitness function, a measure of how well the solution obtained satisfies the requirements of the problem [146]. For the problem of finding a feasible CES variant, we propose the following algorithm for computing the fitness function. For a CES to be able to produce the required amount of product within the given time, it must be operable, i.e., satisfy the compatibility constraint (1.4) and the processing time constraint (1.5). As the ant moves down the decision tree, the number of satisfied compatibility constraints (1.4) for the initial fragments of the CES consisting of stages 1 to i is counted. For example, at tree level i = 1, the initial fragment consists of stage 1; at level i = 2, of stages 1, 2; at level i = 3, of stages 1, 2, 3; at level i = 4, of stages 1, 2, 3, 4; and so on. Operability of the CES is a necessary but not sufficient condition, since the product range must also be produced within the given time, i.e., the processing time constraint (1.5) must be satisfied. Therefore, if constraint (1.5) is satisfied, we add one to the number of satisfied compatibility constraints. Thus, the ideal value of the fitness function (the goal) for a feasible CES variant equals the total number of stages plus one, and we use this value as the constant Q = I + 1. The fitness function Lm makes it possible to assess how close the solution found by a particular ant is to the goal. For the example shown in Fig. 1.2, Q = 5, and the fitness function may range from 1 to 5: an initial fragment consisting of the first stage alone is always operable, so Lm ≥ 1, and Lm = 5 means that all constraints are satisfied and the goal Q has been reached. The fitness values make it possible to compare the paths of different ants. For example, if L1 = 4 for the first ant and L2 = 2 for the second, then the first ant has found a solution closer to the goal. Note that the amount of pheromone deposited, Q/Lm , equals 1 for a feasible solution (Lm = Q) and grows as the fitness value decreases (in this example, Q/L1 = 1.25 and Q/L2 = 2.5). Thus, in the implementation used in all the experiments reported in this book, the deposit is inversely proportional to the fitness value, so the pheromone does not selectively reinforce the paths that are closest to the goal. A deposit that grows with the fitness value (for example, Lm /Q) is a natural alternative; comparing the two rules is a subject for further study. Note that different modifications of ACO may differ in the formula used to update the pheromone concentrations. For example, elitist ants, which have found the best solutions, may be introduced, and only their contributions to the pheromone concentration are taken into account, while all other ants are ignored [147]. In the proposed modification, the pheromone deposited by

57

all ants is taken into account. Over time, the pheromone deposited at previous iterations evaporates, i.e., its concentration decreases (line 10). In the proposed version of the algorithm, the pheromone concentration is updated with evaporation taken into account [148] after the solution of each ant has been evaluated: τrs = τrs · ρ + ∆τrs , where ∆τrs = Q/Lm for the edges of the path of the mth ant and ∆τrs = 0 for the other edges, and ρ ∈ [0, 1] is the evaporation coefficient. Since this update is applied after each ant, evaporation acts M times per iteration of the colony, and the pheromone values reflect mainly the paths of the most recently evaluated ants. The first version of the solution of the problem used the value ρ = 0.9.

5.4

GPU Implementation of the Hybrid Algorithm (ACO+SA)

Figure 5.1 illustrates the five-step implementation of the hybrid algorithm (ACO + SA) on computing systems with GPUs.

Figure 5.1: Structure of the hybrid algorithm (ACO + SA) 1. The CPU reads the input data, initializes the parameters of ACO and SA, transfers all necessary data to the GPU, and then launches the kernel function for ant colony optimization (the ACO kernel). 2. The ACO kernel performs a parallel search for a feasible CES variant as described in Section 5.3. It uses the Multiple Ant Colonies (MAC) method [149]: all colonies run as parallel threads and solve the problem independently of one another. As soon as any colony finds a solution, the computation stops. Not all parallel threads need to find a feasible solution; one is enough. Thus, as the number of parallel threads increases, the probability of finding a feasible solution increases and, accordingly, the search time decreases. 3. The CPU receives the feasible initial solution found, distributes it among the threads, and launches the kernel function for simulated annealing (the SA kernel). 4. The SA kernel searches for the optimal CES variant by simulated annealing, using the feasible CES variant obtained at the previous step as the starting point. It uses Multiple Independent Runs (MIR) [150], in which every parallel thread must find its own feasible CES variant. The threads exchange only the cost minCost of the best solution found, which is used in the acceptance criterion for candidate solutions. The run time is practically independent of the number of threads, but as the number of threads increases, so does the probability that some thread will find the global optimum or a solution close to it. The kernel function finishes when every thread has completed its search. 58

5. The CPU receives the solutions obtained by all threads at the previous step and selects the best of them. This best solution is the final solution of the problem. The pseudocode of the algorithm implementing the proposed approach for the host (CPU) and the device (GPU) is shown in Listings 5.3, 5.4 (ACO), and 5.5 (SA), respectively. Most variable names coincide with those used in the problem statement. Host (CPU) The host (CPU) starts by calling the procedure ReadInputData(), which reads the input data from a file (line 2). The variable inData is a pointer to a structure whose fields are variables and arrays of various dimensions holding all the necessary input data. These data are the same for all threads and do not change while the problem is being solved, i.e., they are constants. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

main () { /* prepare the input data */ ReadInputData ( inData ) ; /* the number of threads is a command - line p a r a m e t e r */ numThreads = argv [1]; /* prepare the working data */ P r e p O p e r a t i o n a l D a t a ( oprData , inData , numThreads ) ; /* t r a n s f e r the n e c e s s a r y data to the device */ cud aMemcpyH toD ( inData , oprData ) ; /* d e t e r m i n e the kernel launch p a r a m e t e r s */ blocksPerGrid = numThreads / M A X _ T H R E A D S _ P E R _ B L O C K + 1; th re ad s Pe rB lo ck = M A X _ T H R E A D S _ P E R _ B L O C K ; /* i n i t i a l i z e the array t h r e a d I d x with ’ -1 ’ */ for ( n = 0; n < numThreads ; n ++) { threadIdx [ n ] = -1; } isFound = 0; /* i n i t i a l i z e the global thread control flag for ACO */ cuda MemcpyHt oD ( isFound , threadIdx ) ; do { /* launch the kernel to find a f e a s i b l e initial s o l u t i o n */ ACO < < < blocksPerGrid , threadsPerBlock > > >() ; c u d a D e v i c e S y n c h r o n i z e () ; /* s y n c h r o n i z e with the device */ cuda MemcpyDt oH ( isFound ) ; } while ( isFound == 0) ; /* copy the f e a s i b l e initial s o l u t i o n from the device to the host */ cuda MemcpyDt oH ( Wfirst , threadIdx ) ; /* find the ID of the thread that has found a f e a s i b l e s o l u t i o n */ for ( n = 0; n < numThreads ; n ++) { if ( threadIdx [ n ] != -1) { tid = threadIdx [ n ]; break ; } } /* the f e a s i b l e s o l u t i o n is the initial a p p r o x i m a t i o n for all threads */ for ( n = 0; n < numThreads ; n ++) { W [ n ] = Wfirst [ tid ]; } /* copy the initial a p p r o x i m a t i o n to the device */ cuda MemcpyHt oD ( W ) ; /* launch the kernel to find the optimal s o l u t i o n */ SA < < < blocksPerGrid , threadsPerBlock > > >() ; c u d a D e v i c e S y n c h r o n i z e () ; /* s y n c h r o n i z e with the device */ /* copy the search results from the device to the host */ cuda MemcpyDt oH (W , Cost , minCost ) ; /* find the best s o l u t i o n among all s o l u t i o n s found */ for ( n = 0; n < numThreads ; n ++) { if ( Cost [ n ] <= minCost ) { Wopt = W [ n ]; minCost = Cost [ n ]; } } }

Listing 5.3: Pseudocode of the host (CPU) algorithm The number of threads used for the computation, numThreads, is a parameter of the program and is passed to it as a command-line argument (line 4). This value is then used when calling the procedure PrepOperationalData(), which prepares the working data (line 6). The variable oprData is a pointer to a structure whose fields are variables and arrays of various dimensions whose values are modified independently by all threads while the problem is being solved. Therefore, one of the array indices is the thread identifier threadID. The procedures 59

ReadInputData() and PrepOperationalData() perform all memory allocations needed for the arrays and initialize the variables, including the pseudorandom number generators of all threads. The data are then transferred to the device for use in subsequent kernel calls (line 8). Next, the kernel launch parameters are determined, i.e., the required number of blocks in the grid, blocksPerGrid (line 10), and the number of threads per block, threadsPerBlock (line 11), are computed. When the required number of blocks is computed (line 10), the result of the division may be non-integer and is automatically truncated to an integer; therefore, one is added to the quotient to round it up. Next, the host initializes the elements of the array threadIdx with ’-1’ (line 13). The length of threadIdx equals the total number of threads numThreads. If a thread with index tid finds a feasible solution, the index of this thread is written to the array: threadIdx[tid] = tid. Knowing the index of the thread that has found a solution, one can then retrieve the feasible solution found from the array of solutions, Wfirst[tid]; this solution is in turn used as the feasible initial solution for simulated annealing. The thread control flag isFound is used to control the threads on the GPU. The flag is initialized to zero (line 14). If isFound==0, a thread continues working. Otherwise, the first feasible solution has been found, isFound!=0, and all threads in the grid stop. The host then transfers the data to the device (line 15) and launches the kernel function ACO() (line 17). The launch is repeated until at least one thread has found a feasible solution, i.e., as long as the value of isFound, copied from the device after each launch, is zero (lines 16–20). In CUDA, a kernel function is launched asynchronously: control returns to the CPU immediately after the kernel has been launched, without waiting for it to finish on the GPU. By calling cudaDeviceSynchronize() (line 18), the host pauses and waits until all threads on the GPU have finished (i.e., synchronizes with the device). After the GPU has finished, the results of the search for a feasible solution are transferred to the host (line 22). Here Wfirst is a two-dimensional array of numThreads vectors of length I each. The elements of each vector hold the selected standard size of the unit for each stage. Having received the solution, the host determines the identifier tid of the thread that has found it by searching for the first element of threadIdx not equal to ’-1’ (lines 24–26). The solution Wfirst[tid] with this identifier is then used as the feasible initial solution for simulated annealing. The host then initializes the initial approximations of all threads with the feasible solution found (lines 28–29). Here W is a two-dimensional array of numThreads vectors of length I each. The elements of each vector hold the selected standard size of the unit for each stage. Next, the host transfers the data to the device (line 31) and launches the kernel function for simulated annealing, SA() (line 33). After synchronization (line 34), the results are copied from the device to the host (line 36). Here Cost is a one-dimensional array of numThreads elements used to store the local costs of the CES variants found by each thread. The host then finds the minimum cost and stores it in the variable minCost (lines 38–41). The optimal CES variant is stored in the array Wopt with I elements, each holding the standard size of the unit for the corresponding stage. Device (GPU) In the proposed parallel implementation, each thread simulates the search of one ant colony (swarm) consisting of M ants [151, 152]. In particular, M = 100 was used in the computational experiments. There is no parallelization over individual ants: each colony searches for its own solution independently. Pseudorandom numbers are generated with NVIDIA’s specialized CUDA Random Number Generation library, cuRAND [153]. Its function curand_uniform() returns a uniformly distributed random value between 0.0 and 1.0, excluding the lower bound (0.0) and including the upper bound (1.0). The cuRAND library is specifically optimized for NVIDIA GPUs and is widely used in implementations of stochastic algorithms on GPUs [154, 155]. 60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

_ _ g l o b a l _ _ void ACO () { /* get the thread i d e n t i f i e r */ threadID = blockDim . x * blockIdx . x + threadIdx . x ; /* if t h r e a d I D is less than the number of threads */ if ( threadID < numThreads ) { /* i n i t i a l i z e the p h e r o m o n e matrix */ for ( i = 1; i <= I ; i ++) { for ( j = 1; j <= J [ i ]; j ++) { tau [ i ][ j ] = curand_u niform () ; } } iterCounter = 0; /* i n i t i a l i z e the local i t e r a t i o n counter */ bestAntId = -1; /* until a s o l u t i o n has been found */ while ( isFound == 0 && iterCounter < m a x I t e r a t i o n N u m b e r ) { /* c o n s t r u c t a path for each ant m of the colony of size M */ for ( m = 1; m <= M ; m ++) { for ( i = 1; i <= I ; i ++) { sum = 0.0; for ( j = 1; j <= J [ i ]; j ++) { eta [ i ][ j ] = ( S ( Want [ m ] , i ) == 0 ? 1:0) + X [ i ][ j ]; sum += pow ( tau [ i ][ j ] , alpha ) * pow ( eta [ i ][ j ] , beta ) ; } r = curand_u niform () ; sump = 0.0; for ( j = 1; j <= J [ i ]; j ++) { p = pow ( tau [ i ][ j ] , alpha ) * pow ( eta [ i ][ j ] , beta ) / sum ; sump += p ; if ( sump > r ) { Want [ m ][ i ] = j ; break ; } } } } /* compute the change in p h e r o m o n e c o n c e n t r a t i o n */ for ( m = 1; m <= M ; m ++) { L [ m ] = NumberS ( Want [ m ] , I ) + ( T ( Want [ m ]) <= Tmax ? 1:0) ; if ( L [ m ] == Q ) { atomicAdd ( isFound , 1) ; bestAntId = m ; } for ( i = 1; i <= I ; i ++) { for ( j = 1; j <= J [ i ]; j ++) { dtau [ i ][ j ] = 0.0; } } for ( i = 1; i <= I ; i ++) { idx = Want [ m ][ i ]; dtau [ i ][ idx ] += Q / L [ m ]; } /* update the p h e r o m o n e c o n c e n t r a t i o n with e v a p o r a t i o n */ for ( i = 1; i <= I ; i ++) { for ( j = 1; j <= J [ i ]; j ++) { tau [ i ][ j ] = tau [ i ][ j ] * rho + dtau [ i ][ j ]; } } } iterCounter ++; } /* i n c r e m e n t the i t e r a t i o n counter */ /* store the first f e a s i b l e s o l u t i o n and the ID of the thread that found it */ if ( bestAntId != -1) { Wfirst [ threadID ] = Want [ bestAntId ]; threadIdx [ threadID ] = threadID ; } } } _ _ d e v i c e _ _ int NumberS (W , level ) { count = 0; for ( i = 1; i <= level ; i ++) { if ( PartS (W , i ) == 0) count ++; } return count ; }

Listing 5.4: Pseudocode of the ant colony optimization (ACO) kernel

The kernel function ACO() starts by determining the thread identifier threadID (line 3). Its value must be less than the required number of threads numThreads (line 5). This check is necessary because the number of blocks in the grid is rounded up when the kernel parameters are computed (see Listing 5.3, line 10). The number of threads actually launched on the GPU may exceed numThreads, whereas all dynamic memory allocations for the arrays in oprData were made for exactly numThreads elements. For all threads whose identifier is less than numThreads, the matrix of pheromone concentrations is initialized with pseudorandom values between 0 and 1 (lines 7–8), and the local iteration counter iterCounter is reset to zero (line 9). The global flag isFound and the local iteration counter iterCounter are used to control the threads: each thread works only while isFound==0 and iterCounter<maxIterationNumber (line 12). The value of the global flag isFound can be changed by any thread that has found a feasible solution, using the atomic function atomicAdd() [156]. Atomicity guarantees that the operation is performed correctly when several threads attempt to execute it simultaneously [109]. The local iteration counter iterCounter protects the algorithm from looping forever if no 61

feasible solution can be found. If no ant of the colony manages to find a feasible solution within maxIterationNumber iterations (which is quite possible for a stochastic algorithm), the thread terminates. After initialization, a path from the first to the last level of the tree is built for each ant m of the colony of size M (i.e., solutions of the problem are constructed) (lines 14–23). Here Want is a local two-dimensional array of M elements, each of which is a vector of length I. The elements of each vector hold the selected main size of the unit for each stage. The heuristic information eta[i][j] is computed as the sum of the operability indicator of the initial fragment of the CES variant consisting of stages 1 to i (line 17) and the main size of the unit X[i][j]. The pheromone values tau[i][j] and the heuristic values eta[i][j] are stored per stage i and unit j, i.e., they are shared by all edges of level i of the search tree that correspond to the same unit. The probability is then computed (line 21), and the node to move to is selected (lines 19–23) as described in Section 5.3. The pheromone concentrations are changed according to the value of the fitness function, which measures the quality of the solution found by the mth ant. The fitness value L[m] (line 26) of each ant m (line 25) equals the number of stages in the solution for which constraint (1.4) is satisfied, plus 1 if constraint (1.5) is satisfied for this solution. The number of stages in the solution for which constraint (1.4) is satisfied is determined by the function NumberS() (lines 42–47). If L[m]==Q (a solution has been found, line 27), the global flag isFound is changed and the index of the ant that has found the solution is stored in bestAntId (line 27); the kernel ACO() then stops at the next iteration of the algorithm. Next, the changes in pheromone concentration dtau[i][j] are computed (lines 28–31). At the next step, the pheromone concentrations tau[i][j] are updated with evaporation taken into account (lines 33–35), after which the iteration counter iterCounter is incremented (line 36). If a feasible solution has been found (line 38), it is stored in Wfirst (line 39), and the thread index threadID is stored in the array threadIdx (line 40). The kernel function for simulated annealing, SA() (Listing 5.5), starts by determining the thread identifier threadID (line 3), in the same way as the ant colony kernel (Listing 5.4). Its value must be less than the required number of threads numThreads (line 5). For all threads whose identifier is less than numThreads, the current local solution Wloc is initialized with the feasible solution W[threadID] found by ACO (line 7), and the temperature t is initialized with the initial value Tinit (line 7). Here Wloc and Wcand are one-dimensional arrays of I elements, each of which stores the standard size of the unit for the corresponding stage. The search continues until the temperature t drops below the final temperature Tfinal (line 8). Before the current temperature t is lowered, the process is repeated Lmax times (line 9). At each iteration, a candidate solution Wcand is constructed in the neighborhood of the current solution Wloc by the procedure Perturb() (line 11). To do so, Perturb() (lines 30–36) randomly chooses a stage stage in the current feasible solution (line 32) and then, for this stage only, randomly chooses a unit from the set J[stage] of units admissible for it (line 34). That is, at each iteration we replace one unit at one stage of the current feasible solution. The candidate solution Wcand is checked for feasibility, i.e., for satisfying the compatibility constraint (1.4) and the processing time constraint (1.5) (line 13). If the candidate is feasible, the difference deltaCost between its cost and the cost minCost of the best solution found so far by all threads is computed (line 15). If deltaCost is negative, the candidate is accepted with probability p = 1; otherwise, it is accepted with probability p = exp(-deltaCost/t), which depends on deltaCost and on the temperature t (line 17). If p exceeds a random number between 0 and 1 (line 19), the candidate Wcand becomes the current solution Wloc (line 21) and is used for the further search. If the cost of the new current solution Wloc is lower than the cost minCost of the current

62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

_ _ g l o b a l _ _ void SA () { /* get the thread i d e n t i f i e r */ threadID = blockDim . x * blockIdx . x + threadIdx . x ; /* if t h r e a d I D is less than the number of threads */ if ( threadID < numThreads ) { /* i n i t i a l i z e the current s o l u t i o n with the f e a s i b l e initial s o l u t i o n */ Wloc = W [ threadID ]; t = Tinit ; while ( t > Tfinal ) { /* until the final t e m p e r a t u r e is reached */ for ( l = 0; l < Lmax ; l ++) { /* c o n s t r u c t a new c a n d i d a t e s o l u t i o n */ Wcand = Perturb ( Wloc ) ; /* check the c o m p a t i b i l i t y and p r o c e s s i n g time c o n s t r a i n t s */ if ( S ( Wcand ) == 0 && T ( Wcand ) <= Tmax ) { /* compute the change in the o b j e c t i v e f u n c t i o n */ deltaCost = DefineCost ( Wcand ) - minCost ; /* compute the p r o b a b i l i t y */ p = deltaCost < 0 ? 1.0 : exp ( - deltaCost / t ) ; /* if the p r o b a b i l i t y exceeds a random number */ if ( p > cu rand_uni form () ) { /* accept the new s o l u t i o n */ Wloc = Wcand ; /* if the s o l u t i o n found is better */ if ( DefineCost ( Wloc ) < minCost ) { W [ threadID ] = Wloc ; Cost [ threadID ] = DefineCost ( Wloc ) ; atomicMinDbl ( minCost , Cost [ threadID ]) ; }}}} t *= sigma ; } } /* update the t e m p e r a t u r e */ } _ _ d e v i c e _ _ int * Perturb ( W ) { Wn = W ; stage = ( int ) curand (1 , I ) ; /* choose a random stage */ /* choose a random s t a n d a r d unit size for this stage */ Wn [ stage ] = ( int ) curand (1 , J [ stage ]) ; return Wn ; }

Listing 5.5: Pseudocode of the simulated annealing (SA) kernel

global optimal solution (line 23), the new solution replaces the current global optimal solution W (line 24), and its cost is stored in the array Cost (line 25) and in the global variable minCost using the atomic function atomicMinDbl() (line 26). The temperature t is then lowered (line 28). If no feasible candidate solution is obtained during the whole run, as the temperature decreases from its maximum to its minimum value, the result is the feasible initial solution used as the starting point (line 7).

5.5

Results of Computational Experiments

For the computational experiments, a program was developed in C using CUDA [109]. The experiments were carried out on a computing system with an Intel Xeon processor (E5-1620 v2, 4 cores with Hyper-Threading, 10 MB Smart Cache, 3.7 GHz) and 16 GB of RAM, equipped with an NVIDIA Tesla K20c GPU. GPU specifications: 13 streaming multiprocessors (SMs) with 192 CUDA cores each (2496 CUDA cores in total), 5 GB of global memory, and 48 KB of shared memory per SM. The software comprised Ubuntu 16.04.1, NVIDIA driver 367.57, CUDA 8.0, and the GNU C++ compiler 5.4.0. The test case was the CES equipment selection problem described in Section 1.3 (see Fig. 1.4), with 2 to 12 standard equipment sizes available for each of the 16 processing stages (216 to 1216 CES variants in total). In Chapters 2–4, the optimal CES equipment selection problem was solved by branch-andbound, which yields globally optimal solutions [54, 94, 95]. Here, on the same test case and the same test platform, we compare the solutions obtained by the hybrid metaheuristic algorithm 63

(ACO + SA) with those obtained by branch-and-bound. Since both ant colony optimization and simulated annealing are stochastic (probabilistic) algorithms, the results for the same problem may differ even for identical input data. Therefore, for each example, the algorithm was run 100 times and the results were averaged. The first experiment studied the effect of the number of threads on the program’s performance (see Fig. 5.2). 450

400

Run time, s

350

300 250

ACO run time SA run time

200 150

100

2500

2400

2300

2200

2100

1900

Number of threads

2000

1800

1700

1500

1600

1400

1300

1200

1100

1000

900

800

700

500

600

400

300

200

0

100

50

Figure 5.2: Run time of the parallel program versus the number of threads The experiments were carried out for the test CES with 10 standard equipment sizes available at each stage (1016 possible CES variants in total). The program was run with 100 to 2500 threads in steps of 100. The run time of the algorithm decreases as the number of threads increases. With 100 threads, the ant colony algorithm (ACO) accounts for 91% of the total run time; this share then decreases and, beyond 2100 threads, remains practically constant at about 11% of the total run time. Figure 5.3 shows the relative deviation of the objective function value (the cost of the CES variant) versus the number of threads used. The relative deviation was computed as deviation =

|observed − expected| · 100%, expected

where the observed value is the cost of the CES variant obtained by the metaheuristic algorithm and the expected value is the globally optimal cost of the CES variant obtained by branchand-bound. With 100 threads, the relative deviation is about 0.38%; starting from about 2000 threads, it decreases practically to 0%. A number of tests for problems of different sizes were carried out to measure the objective function values and run times of the proposed hybrid metaheuristic algorithm (sequential CPU and parallel GPU versions) and of the sequential branch-and-bound algorithm (B&B). The results for the different approaches are shown in Figs. 5.4 and 5.5 (the vertical axis is logarithmic). Figure 5.4 shows the relative deviation of the objective function value (the cost of the CES variant) versus the problem size. The program was run on the GPU with 2500 parallel threads. 64

0.40%

0.30% 0.25% 0.20% 0.15% 0.10%

2500

2300

2400

2100

2200

2000

1800

Number of threads

1900

1700

1600

1500

1300

1400

1100

1200

1000

900

800

700

600

500

300

200

0.00%

400

0.05%

100

Relative deviation, %

0.35%

Figure 5.3: Relative deviation of the objective function value versus the number of threads The parallel version “(ACO+SA) on GPU” obtains near-optimal solutions (for problems of size 216 –916 , the relative deviation is 0%; for problems of size 1016 –1216 , it is less than 0.01%). The relative deviation of the sequential version “(ACO+SA) on CPU” is below 1% for a small problem (216 ), but then increases to 14% for the largest problem (1216 ). This is because a single run of simulated annealing has a low probability of finding a good solution for large problems. Figure 5.5 shows the run time of the program versus the problem size (the vertical axis is logarithmic). The run time of branch-and-bound grows exponentially with the problem size. For the smallest problem (216 possible CES variants in total), the run time is less than 1 second. For a large problem (1216 ), however, it is about 134 hours. The run time of the sequential hybrid algorithm is less than 0.1 second for a small problem (216 ) and increases to 581 seconds for the largest problem (1216 ). The run time of the parallel hybrid algorithm (ACO + SA) increases smoothly from 5 to 227 seconds. The parallel hybrid algorithm is slower than the sequential one for problems with 216 to 916 variants. For the largest problems, with 1016 to 1216 variants, the parallel implementation is approximately 2.6–2.9 times faster than the sequential one. At first glance, a speedup of 2.6–2.9 over the sequential version is modest. However, the quality of the solutions obtained by the parallel implementation “(ACO+SA) on GPU” is higher than that of the sequential “(ACO+SA) on CPU”: the relative deviation of the solution obtained from the globally optimal solution found by branch-and-bound is less than 0.01%, compared with 9–14% for the sequential program. This is achieved because the parallel version performs many simultaneous independent searches for the optimal solution, and the more parallel threads there are, the more likely it is that, in practically the same run time, some thread will find the globally optimal solution. Running the parallel algorithm with 2500 threads is equivalent to running the sequential algorithm 2500 times and choosing the best of the solutions found. The experiments confirm that the parallel implementation of the hybrid metaheuristic algorithm (ACO+SA) yields high-quality solutions for the test problem that coincide with the globally optimal values obtained by the deterministic branch-and-bound algorithm, but much 65

100.0%

(ACO + SA) on CPU

Relative deviation, %

10.0%

(ACO + SA) on GPU

1.0% 0.1% 0.01%

0.001%

0.0001%

216

316

416

516

616

716

816

916

1016

Number of CES variants (problem size)

1116

1216

Figure 5.4: Relative deviation of the objective function value versus the problem size faster: in a few tens of seconds instead of a few tens of hours.

66

1.0E6 1.0E5

1.0E4

B&B on CPU (ACO+SA) on CPU (ACO+SA) on GPU

Run time, s

1.0E3

100.0

10.0 1.0

0.1

0.01

0.001 216

316

416

516

616

716

816

916

Number of CES variants (problem size)

1016

Figure 5.5: Run time of the program versus the problem size

67

1116

1216

Chapter 6

Tuning the Parameters of the Ant Colony Algorithm The hybrid algorithm presented in Chapter 5 works in two stages: (1) finding a feasible initial solution (a fully operable CES variant that produces the given product range within the given time) by ant colony optimization (ACO), and (2) improving this solution by simulated annealing (SA). The first stage is thus the basis for the second, and without it the hybrid algorithm cannot work. In Chapter 5, a feasible initial solution was successfully found by ant colony optimization for problems with up to 1216 (≈ 1017 ) possible CES variants [157, 158]. However, when the problem size exceeds 1116 , the performance of the algorithm deteriorates considerably despite the use of high-performance GPUs. Problems larger than 1216 were not solved, which motivated us to turn to tuning the parameters of the algorithm. The chapter is organized as follows. Section 6.1 reviews existing approaches to tuning the parameters of metaheuristic algorithms. Section 6.2 discusses the effect of the parameters α, β, and ρ on the ant colony algorithm. Section 6.3 proposes a method for tuning these parameters based on statistical analysis of the results of multiple runs of the parallel program, together with three ways of applying it: random, constant, and multi-constant. Section 6.4 reports the results of computational experiments and compares the approaches in terms of the time needed to find a feasible solution and the time spent on the tuning study. Section 6.5 presents the results of running the program on GPUs of three generations (Kepler, Pascal, Volta) and on compute nodes with several GPUs.

6.1

Approaches to Tuning the Parameters of Metaheuristic Algorithms

It is well known that the efficiency of metaheuristic algorithms depends to a large extent on their parameters [26, 159]. Identifying these parameters (so-called tuning of metaheuristics) is a nontrivial problem to which a considerable body of research has been devoted. As a rule, there are no universally optimal parameter values for metaheuristic algorithms [8]. There are two strategies for setting the parameters of an algorithm: offline parameter initialization and online parameter initialization [8, 160]. With offline initialization, the parameter values are searched for before the algorithm is actually applied. The best values are determined once at the design stage and are not changed thereafter when the algorithm is used [161]. Traditionally, offline tuning has mainly been performed in a trial-and-error fashion [148]. The drawback of this approach is that trial-and-error tuning depends heavily on the intuition and experience of the algorithm designer, which are usually difficult or impossible to document and, therefore, to reproduce. Recently, much attention has been paid to search-based offline

68

tuning methods, which allow the parameters of algorithms to be tuned (semi-)automatically with minimal human involvement. Nevertheless, these methods yield specific parameter values before the algorithm is applied, and these values do not change while a particular problem is being solved. An alternative to offline tuning is online tuning [162]. Typically, online tuning makes it possible to modify the parameters of the algorithm while a particular problem is being solved. The potential advantage of online tuning is that algorithms can adapt to the characteristics of a particular problem and solve it more efficiently. When the problems to be solved are relatively heterogeneous, parameter values found in advance, which generally lead to good performance on average, may lead to considerably worse results in some individual cases. Adapting the parameters during the run may also be useful for achieving the best performance depending on the stage of the search. The operation of an algorithm can often be divided into explorative and exploitative search phases, and the optimal parameter values in these phases may differ considerably. Furthermore, if algorithms are applied in situations very different from those for which they were designed or tuned offline, adapting the parameters during the run increases the robustness of the algorithm. There are various ways of changing the parameters of an algorithm during its run [161]. One of the simplest is to specify a rule for changing the parameters before the algorithm is actually launched. In this approach, the static parameters are replaced offline by (deterministic or stochastic) functions of the computation time or of the number of iterations. Some authors use the term deterministic parameter control [163]. Others, however, consider this term inaccurate and propose the term pre-scheduled parameter variation instead, since the adjective “deterministic” does not quite characterize this method of parameter control, given that the parameters may be chosen at random [164]. Although pre-scheduled parameter variation is an online tuning method, it does not eliminate the problem of offline tuning, because in this case, too, the parameters must be set appropriately before the algorithm is launched. An alternative approach is to use adaptive parameter settings, in which the scheme for changing the parameters is defined as a function of some statistics of the algorithm’s behavior [165]. The indicators used for online adaptation can be roughly grouped according to whether they are based on absolute or relative information. In the first case, the adaptation strategy monitors the occurrence of certain events during the run, for example, reaching some fixed threshold distance between the solutions found. Exceeding this threshold triggers a set of predefined rules for changing the parameters of the algorithm. In the second case, the adaptation strategy takes into account the relative difference between the performance achieved with different parameter values and adapts the parameter values towards the most successful ones. Parameter adaptation requires certain design decisions that are not directly related to the implementation of the algorithm itself. In particular, the expressions describing the parameter updates must be defined a priori, and one must be confident that the mechanisms used are very robust with respect to their particular implementation. Another object of study may be tuning parameters that are modified by the algorithm itself during its run. For example, the tuning parameters of the search algorithm are added to the search space of the main problem being solved, and the optimization problem is then solved in the new, extended search space. This approach is called self-adaptation [166–168]. Some authors extend the notion of self-adaptation to search-based adaptation, in which the parameters of the algorithm are adapted using approaches different from the one used to solve the main problem [159, 169]. Moreover, all existing approaches to studying the relationship between parameter settings, problem characteristics, and algorithm performance can be divided into three categories: (1) analytical approaches, (2) automated approaches, and (3) empirical approaches [170, 171]. The analytical approach to understanding algorithm parameters and recommending their values is based on mathematical proofs. The automated approach consists in applying another algorithm or heuristic to tune the parameters of the main algorithm. If such an automated

69

approach uses the algorithm or heuristic in an introspective sense (for self-analysis during the run), the term automated self-tuning is often used. The empirical approach collects experimental data on the behavior of the algorithm and analyzes them to draw conclusions by trial and error or by the One-Factor-At-a-Time (OFAT) method, and to give recommendations on choosing parameter values. Design of Experiments (DOE) methods are typically used for this purpose [172, 173]. Among automated offline tuning methods, the most widely used are so-called algorithm configurators, general-purpose software tools that treat the algorithm being tuned as a black box and search for its best configuration on a training set of problem instances within a given computational budget [27, 160]. Configurators differ in their search strategy. Racing methods, in particular F-Race [174] and its iterated version irace [175], generate a set of candidate configurations, evaluate them sequentially on more and more problem instances, and use statistical tests (e.g., the Friedman test) to discard clearly inferior candidates early. The distributions from which the candidates of the next iteration are generated gradually concentrate around the best (“elite”) configurations. ParamILS [176] performs an iterated local search in the space of discretized parameter values, while SMAC [177] builds a regression model of algorithm performance as a function of the parameter values and instance features and uses it to select the most promising configurations. What these methods have in common is that they are designed for expensive sequential runs of the algorithm being tuned: their main goal is to reduce the number of runs needed to find a good configuration.

6.2

Parameters of the Ant Colony Algorithm

Before the ant colony algorithm is started, a number of parameters that affect its operation must be set [178]. These parameters have an extremely strong influence on the quality of the solutions obtained and on the run time of the algorithm [179]. The parameters α (the information elicitation factor) and β (the expected heuristic factor) control the balance between the influence on the search of the pheromone concentration (which reflects the history of successful solutions) and of the heuristic information (which expresses the desirability of a particular choice based on the features of the problem), respectively. For α = 0, the pheromone information is not used, and the previous experience of the ant colony is ignored. This reduces the convergence rate of ACO but increases its capacity for global search. If α is large, an ant will predominantly choose the same path as the previous ants [179]. The convergence rate of ACO then increases, but so does the probability of getting trapped in a local optimum, while the capacity for global search decreases. Conversely, if α is too small, the convergence of ACO slows down, but the global search capability of the algorithm increases. If β = 0, the heuristic information is not taken into account, and an ant chooses its path solely on the basis of the pheromone concentration, which leads to suboptimal solutions. If β is large, the heuristic information dominates the choice of path, and the algorithm approaches a greedy algorithm [180, 181]. For example, in the travelling salesman problem, the city nearest to the current one will be chosen [182]; in the search for a feasible initial CES variant, the unit of the largest size will be chosen. The pheromone evaporation rate ρ ∈ (0, 1] controls the degree to which the pheromone concentration decreases. In the literature, the expression (1 − ρ) is also used as the factor by which the current pheromone concentration is multiplied upon evaporation [142]. This can be interpreted as follows: as the evaporation rate ρ increases, the factor (1 − ρ) decreases, and the pheromone evaporates faster: a value close to one means strong evaporation, and a value close to zero means weak evaporation. For example, for ρ = 0.9 (strong evaporation), the factor is 1 − ρ = 0.1, and multiplying by 0.1 reduces the concentration tenfold. Conversely, for ρ = 0.1 (weak evaporation), multiplying by 1 − ρ = 0.9 reduces the concentration by 10%. In our implementation, the value ρ itself is used as the factor, so the logic is reversed: a value of ρ close to one means weak evaporation, and a value close to zero means strong evaporation [148]. 70

As noted in Section 4.3, differences in the control flow of individual threads in a grid degrade overall performance because of warp divergence. Different numbers of ants in the colonies would lead to different numbers of loop iterations in the GPU threads, which would degrade the computational performance of the GPU. For this reason, the number of ants is not tuned in this book; it is the same for all colonies and equal to 100.

6.3

Proposed Tuning Method

The goal of tuning is to find values of the parameters α, β, and ρ for which ant colony optimization finds a feasible initial solution as fast as possible. Here a constraint satisfaction problem is solved, i.e., the goal is to find a feasible solution rather than a solution that is optimal in terms of capital cost. The parameters are tuned by statistical analysis of the results of multiple experimental runs of the program. We analyze how frequently solutions are found for particular intervals of the parameters α, β, and ρ. Checking how often successful results occur within a range of parameter values does not require any specific and non-obvious information, such as the shape of membership functions in fuzzy logic [183, 184] or automated tuning methods [185, 186]. Owing to the parallelization of the algorithm and the use of modern GPUs, a large number of computational experiments can be performed to accumulate statistics in acceptable time. The idea of the proposed method is close to iterated racing [175]: in both cases, the search region for the parameter values is successively narrowed around the most successful values. However, there are several differences. First, the evaluation of configurations is built into the parallel algorithm itself. Each GPU thread is an independent ant colony with its own triple of values α, β, and ρ; therefore, a single run of the program tests thousands of configurations simultaneously, and no external control software is needed for tuning. Second, the success criterion for a configuration is not averaged performance but the very fact that a feasible solution has been found. For a thread that has found a feasible solution, its triple of values is stored, and the intervals are narrowed according to the frequency distribution of such triples. This matches the constraint satisfaction problem being solved, in which the result of each run is effectively binary: a feasible solution is either found or not. Third, the multi-constant approach produces not a single “best” configuration but a parallel portfolio of configurations [187] executed simultaneously by different threads, which is natural for the massively parallel GPU architecture. The price paid for the simplicity of the method is the absence of statistical tests for comparing configurations, and the fact that tuning was performed on a single family of problems (one CES with different problem sizes), whereas configurators are usually trained on many heterogeneous instances. The applicability of the parameter values found to other CES requires further verification. We propose three variants of tuning the algorithm: the random approach, the constant approach, and the multi-constant approach [188, 189]. Figure 6.1 shows the general approach to tuning the parameters of the algorithm: 1. The host (CPU) reads the input data and initializes the parameters on the device (GPU). 2. The device initializes the tuning parameters in one of the ways shown in Listing 6.1. 3. The host launches the kernel function ACO() shown in Listing 5.4. 4. The kernel function ACO() searches for a feasible solution. 5. The host receives the search results for further processing; if no feasible solution has been found, step 2 or step 3 is repeated, depending on the tuning approach adopted. Otherwise, the feasible solution found is stored. 71

Figure 6.1: General method for tuning the ACO parameters 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

_ _ g l o b a l _ _ void SET ( TUN_APP ) { ... /* get the thread i d e n t i f i e r */ threadID = blockDim . x * blockIdx . x + threadIdx . x ; switch ( TUN_APP ) { case RANDOM : /* p a r a m e t e r s are i n i t i a l i z e d with random values */ alpha [ threadID ] = rand ( alpha_a , alpha_b ) ; beta [ threadID ] = rand ( beta_a , beta_b ) ; rho [ threadID ] = rand ( rho_a , rho_b ) ; break ; case CONSTANT : /* the same p a r a m e t e r values for all threads */ alpha [ threadID ] = c onstant_ alpha ; beta [ threadID ] = constant_beta ; rho [ threadID ] = constant_rho ; break ; case MULTI_C ONSTANT : /* d i f f e r e n t p a r a m e t e r values for all threads */ alpha [ threadID ] = c onstant_ alpha [ threadID ]; beta [ threadID ] = constant_beta [ threadID ]; rho [ threadID ] = constant_rho [ threadID ]; break ; } ... }

Listing 6.1: Ways of setting the algorithm parameters

Random approach The parameters α, β, and ρ are initialized with uniformly distributed random values from the given intervals [αa , αb ], [βa , βb ], and [ρa , ρb ], respectively (see Listing 6.1, lines 6–8). The interval bounds are taken according to the recommendations of [146]. To generate random numbers in a given interval (a, b] other than (0.0, 1.0], the expression rand(a,b) = curand_uniform() * (b-a) + a is used; because of the specifics of the function curand_uniform(), the lower bound is excluded and the upper bound is included. By running this algorithm many times, we obtain a set of triples of the ACO parameters α, β, and ρ for which feasible solutions were found. The stored parameter values are then processed statistically. If no feasible solution has been found at step 4 of the method after maxIterationNumber iterations, the method returns to step 2 (see Fig. 6.1). Constant approach This approach differs from the random one in that the values of α, β, and ρ are the same for all threads; they are initialized once before the algorithm is started and do not change during its execution (see Listing 6.1, lines 11–13). These initial values are obtained by statistical analysis 72

of the set of parameter triples produced by the random approach, although they could also be taken, for example, from recommendations in the literature [146]. If no feasible solution has been found at step 4 of the method after maxIterationNumber iterations, the method returns to step 3 (see Fig. 6.1). Multi-constant approach This approach differs from the constant one in that the values of α, β, and ρ differ from thread to thread; they are initialized once before the algorithm is started and do not change during its execution (see Listing 6.1, lines 16–18). These values are taken from the set of stored triples found with the random approach. In this case, the random approach in effect trains the parallel ACO algorithm, since the random approach stores only those triples of α, β, and ρ for which a solution was found. If no feasible solution has been found at step 4 of the method after maxIterationNumber iterations, the method returns to step 3 (see Fig. 6.1).

6.4

Results of Computational Experiments

For the computational experiments, a program was developed in C using CUDA [109]. The experiments were carried out on a computing system with an Intel Xeon Gold 5118 processor (12 cores with Hyper-Threading, 2.3 GHz) and 192 GB of RAM, equipped with an NVIDIA Tesla V100-SXM2 GPU. GPU specifications: 80 streaming multiprocessors (SMs) with 64 CUDA cores each (5120 CUDA cores in total), maximum clock rate 1.53 GHz, and 16 GB of global memory. The software comprised CentOS Linux 7.5.1804, NVIDIA driver 410.72, CUDA 10.0, and the GNU C++ compiler 6.4.0. The test case was the CES equipment selection problem described in Section 1.3 (see Fig. 1.4), with 11 to 20 standard equipment sizes available for each of the 16 processing stages (1116 to 2016 CES variants in total). Since both ant colony optimization and simulated annealing are stochastic (probabilistic) algorithms, the results for the same problem may differ even for identical input data. Therefore, for each problem size from 1116 to 2016 (10 series in total), the algorithm was run 100 times; the time needed to find a feasible solution was measured, the average time was computed, and the values of α, β, and ρ were stored. The program was run with 5120 threads, which equals the number of CUDA cores of the GPU used. Random approach The initial search bounds for the first iteration were chosen following the recommendations of [146]: α ∈ (0, 2], β ∈ (0, 2], and ρ ∈ (0, 1]. After a complete series of experiments, 10 · 100 = 1000 triples of the parameter values α, β, ρ were obtained for problems of various sizes. The total GPU time spent on the first series of experiments with the random approach was 51 054 s ≈ 14 h. Figure 6.2 shows the run time of the program versus the problem size. Note that the program does not enumerate all possible CES variants; the search stops at the first feasible solution. Figure 6.3 shows the frequency histogram of the feasible solutions found. The histogram is constructed as follows. The whole interval (0.0, 2.0] of the algorithm parameters is divided into equal subintervals (bins) of width 0.1: (0.0, 0.1], (0.1, 0.2], (0.2, 0.3], . . . , (1.9, 2.0]. These bins are plotted along the horizontal axis, and above each of them a bar is drawn for each of the algorithm parameters α, β, and ρ. The height of a bar is proportional to the number of feasible solutions found with the corresponding parameter value falling into this bin. Analysis of the frequency distribution showed that ≈ 90% of the solutions found were obtained with the ACO parameters in the ranges α ∈ (0.0, 1.2], β ∈ (0.1, 1.2], ρ ∈ (0.2, 1.0]. The idea 73

100 90

α ∈ (0.0, 2.0], β ∈ (0.0, 2.0], ρ ∈ (0.0, 1.0]

80

70

Run time, s

60 50

40 30

20

10

0 1116

1216

1316

1416

1516

1616

1716

1816

Number of CES variants (problem size)

1916

2016

Figure 6.2: Run time of the parallel program versus the problem size at the first tuning iteration of the proposed method for tuning the ACO parameters is to move step by step towards those ranges of the parameters α, β, and ρ in which feasible solutions occur more often (i.e., are found with higher frequency). The search for 1000 new solutions was then repeated with the adjusted parameter ranges, and the frequency of the solutions found was again analyzed for each parameter. This process (finding 1000 solutions, analyzing the frequencies, and adjusting the parameter intervals to those containing 90% of the solutions) was repeated 7 times. The total GPU time after seven tuning iterations was 86 509 s ≈ 24 h. The change in the average run time of the algorithm over the seven iterations for problems of size 1116 –2016 is shown in Fig. 6.4. Figure 6.5 shows the frequency distribution of the feasible solutions found at the last, seventh, iteration (parameter ranges α ∈ (0.1, 0.3], β ∈ (0.3, 0.7], ρ ∈ (0.8, 1.0]). After the last, seventh, iteration, the parameter intervals contracted to the point (α = 0.2, β = 0.5, ρ = 0.9). Thus, the statistical analysis of the random approach yielded parameter values that were subsequently used in the study of the constant approach. Constant approach In the constant approach, every thread running on the GPU uses the same tuple of constants α, β, ρ. Figure 6.6 shows the run time of the program versus the problem size for the constant approach. Table 6.1 gives the average run times of the algorithm over the whole range of problem sizes (from 1116 to 2016 ). The first set of constants, α = 0.4, β = 0.6, ρ = 0.9, was chosen empirically and used in the earlier experiments presented in Section 5.5 [157, 158]. The performance of the algorithm with the new parameter values is compared below with its performance with these values. The second set of constants, α = 0.31, β = 0.54, ρ = 0.65, was taken as the mean values of the parameters obtained after the first iteration of the random approach, which produced 74

450

α ∈ (0.0, 2.0] β ∈ (0.0, 2.0] ρ ∈ (0.0, 1.0]

Number of solutions found

400 350

300 250

200 150

100 50

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

1.1

Interval

1.2 1.3 1.4 1.5

1.6

1.7

1.8 1.9 2.0

Figure 6.3: Frequency distribution of the feasible solutions found after the first iteration Table 6.1: Average run times of the algorithm α

β

ρ

Run time, s

0.4 0.31 0.1 0.2

0.6 0.54 0.3 0.5

0.9 0.65 0.9 0.9

38.1 54.2 4.9 1.5

the first 1000 triples of α, β, and ρ. The intuitive assumption was that averaging known good parameter values would yield good candidate parameters. However, using the averaged parameter values increased the average run time by a factor of ≈ 1.4 compared with the parameters used in Chapter 5 (from ≈ 38.1 s for α = 0.4, β = 0.6, ρ = 0.9 to ≈ 54 s for α = 0.31, β = 0.54, ρ = 0.65). The third set of constants, α = 0.1, β = 0.3, ρ = 0.9, was chosen on the basis of the frequency analysis after the first iteration of the random approach (see Fig. 6.3) as the parameter values approximately corresponding to the highest frequency of feasible solutions. This reduced the average run time of the algorithm by a factor of ≈ 7.7 (from ≈ 38.1 s for α = 0.4, β = 0.6, ρ = 0.9 to ≈ 4.9 s for α = 0.1, β = 0.3, ρ = 0.9). After the seventh iteration, the parameter intervals contracted to the point α = 0.2, β = 0.5, ρ = 0.9, which was used as the fourth set of constants. This reduced the average run time of the algorithm by a factor of ≈ 25 (from ≈ 38.1 s for α = 0.4, β = 0.6, ρ = 0.9 to ≈ 1.5 s for α = 0.2, β = 0.5, ρ = 0.9). Since the constants were determined from the data obtained after seven iterations of the random approach, the tuning time can be considered the same as for the random approach, 86 509 s ≈ 24 h.

75

100 90

80

70

Run time, s

60

Iteration 1: α ∈ (0.0, 2.0], β ∈ (0.0, 2.0], ρ ∈ (0.0, 1.0] Iteration 2: α ∈ (0.0, 1.2], β ∈ (0.1, 1.2], ρ ∈ (0.2, 1.0] Iteration 3: α ∈ (0.0, 0.8], β ∈ (0.1, 1.1], ρ ∈ (0.3, 1.0] Iteration 4: α ∈ (0.0, 0.6], β ∈ (0.2, 1.1], ρ ∈ (0.4, 1.0] Iteration 5: α ∈ (0.0, 0.5], β ∈ (0.2, 0.9], ρ ∈ (0.5, 1.0] Iteration 6: α ∈ (0.0, 0.4], β ∈ (0.2, 0.8], ρ ∈ (0.6, 1.0] Iteration 7: α ∈ (0.1, 0.3], β ∈ (0.3, 0.7], ρ ∈ (0.8, 1.0]

50

40 30

20

10

0 1116

1216

1316

1416

1516

1616

1716

1816

Number of CES variants (problem size)

1916

2016

Figure 6.4: Run time of the parallel program versus the problem size over seven tuning iterations Multi-constant approach In the multi-constant approach, each thread running on the GPU uses its own tuple of constants α, β, ρ. For this purpose, we used the set of constants obtained in the first series of experiments with the random approach, for α ∈ (0.0, 2.0], β ∈ (0.0, 2.0], ρ ∈ (0.0, 1.0]. After 100 experiments for the problem size 1116 , we have a set of 100 parameter triples. After the next 100 experiments for the problem size 1216 , we have 200 triples (100 for 1116 + 100 for 1216 = 200); after the next 100 experiments for the problem size 1316 , 300 triples (100 for 1116 + 100 for 1216 + 100 for 1316 = 300); and so on, up to 1 000 triples after the experiments for the problem size 2016 . When the program is run on the GPU (in this case, 5 120 threads were used on the Tesla V100, corresponding to the total number of CUDA cores of this accelerator), the algorithm parameters are initialized cyclically, as described below. Suppose there are 100 unique parameter triples. These 100 unique triples must be used to initialize 5 120 threads; therefore, the parameter values repeat every 100 threads. For 200 triples, they repeat every 200 threads, and so on; for 1 000 unique triples, the constants repeat every 1 000 threads. The run time of the program versus the problem size is shown in Fig. 6.7. In general, the time to solve the problem does not exceed ≈ 7 s. The worst run time (especially for the largest problem, 2016 ) corresponds to 100 unique triples, and the best one to 1000 triples. Starting from 200 unique triples, the average run time of the algorithm differs from that for the maximum of 1 000 triples by a factor of only 1.13 (3.23 s versus 2.86 s). At the same time, the time needed to find 1000 unique triples for problem sizes from 1116 to 16 20 equals the time of the first iteration of the random approach, 51 054 s ≈ 14 h, while the time needed to find a set of 200 triples for problem sizes from 1116 to 1216 is 3 373 s ≈ 56 min, i.e., 15.1 times less.

76

Number of solutions found

600

α ∈ (0.1, 0.3] β ∈ (0.3, 0.7] ρ ∈ (0.8, 1.0]

500 400 300 200 100 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

1.1

Interval

1.2 1.3 1.4 1.5

1.6

1.7

1.8 1.9 2.0

Figure 6.5: Frequency distribution of the feasible solutions found after the seventh iteration Comparison of the approaches Figure 6.8 shows the best results obtained with each approach. For the multi-constant approach, the variant with 200 triples is also shown, since it gave acceptable results with the minimum time spent on finding the parameter values. The fastest algorithm, tuned with the constant approach, α = 0.2, β = 0.5, ρ = 0.9 (average run time ≈ 1.5 s), is ≈ 2.1 times faster than the slowest one, tuned with the multi-constant approach using 200 unique triples (average run time ≈ 3.2 s). Table 6.2: Tuning time for the different approaches Approach

Tuning time

Random Constant Multi-constant (1 000) Multi-constant (200)

86 509 s (≈ 24 h) 86 509 s (≈ 24 h) 51 054 s (≈ 14 h) 3 373 s (≈ 56 min)

However, the time spent on tuning (see Table 6.2) is 86 509 s ≈ 24 h in the first case and 3 373 s ≈ 56 min in the second, i.e., ≈ 26 times less.

6.5

Experiments on GPUs of Different Architectures and on Multiple GPUs

The results presented in Section 6.4 were obtained on a single Tesla V100 GPU. To assess to what extent the conclusions depend on the hardware used, and to explore ways of further reducing the time needed to find a feasible solution, the program was additionally studied on GPUs of three generations and in configurations with several GPUs within a single compute node. In all experiments, as in Section 6.4, the problem of finding a feasible CES variant was solved for the test case (Section 1.3) with problem sizes from 1116 to 2016 variants; for each problem size, the program was run 100 times, and the run time of the ant colony algorithm until the first 77

120

α = 0.4, β = 0.6, ρ = 0.9 α = 0.31, β = 0.54, ρ = 0.65 α = 0.1, β = 0.3, ρ = 0.9 α = 0.2, β = 0.5, ρ = 0.9

100

Run time, s

80

60

40 20 0 1116

1216

1316

1416

1516

1616

1716

1816

Number of CES variants (problem size)

1916

2016

Figure 6.6: Run time of the program for problems of various sizes, constant approach feasible solution was found was measured. The number of threads launched on one GPU was always equal to the number of CUDA cores of the corresponding accelerator.

6.5.1

Organization of Computations on Multiple GPUs

In the parallel implementation of the ant colony algorithm (Section 5.4), each GPU thread simulates an independent ant colony, and the search stops as soon as one of the colonies finds a feasible solution. The colonies do not exchange information; therefore, extending this scheme to several accelerators requires no changes to the kernel function ACO(): it suffices to launch it simultaneously on all GPUs of the node and to make sure that the search stops on all devices once a solution has been found on one of them. Thus, the total number of simultaneously running colonies increases in proportion to the number of accelerators. The pseudocode of the corresponding part of the host program is shown in Listing 6.2. In the preparatory phase (lines 2–6), the input and working data are copied to each device, after which the pseudorandom number generators and the parameters of the ant colony algorithm are initialized. The generator seed is multiplied by the device number (line 5) so that the random number sequences, and hence the behavior of the colonies, differ between GPUs. The parameters α, β, and ρ are set in one of the ways described in Section 6.3 (line 6; see Listing 6.1). With the random approach, each device generates its own parameter values from the given intervals. With the multi-constant approach, each device is assigned its own contiguous block of triples from a common file, so that all threads of all devices receive different triples. The host then launches the kernel function ACO() on all devices (lines 8–11). The launches are asynchronous, each in its own CUDA stream, and after each launch an event (CUDA event) marking the completion of the kernel is placed into the same stream. The host then polls the state of the events in a loop without blocking on any single device (lines 13–22). As soon as the kernel on some device has finished, the host reads the flag isFound from it. If a solution has been found (line 18), the flag value is copied asynchronously to all devices through separate streams. Since the continuation condition of the main loop of the kernel includes a check of this flag (see Listing 5.4, line 12), the colonies on the other devices stop at their next iteration. After all 78

9

100 triples (α, β, ρ) for problem size 1116 200 triples (α, β, ρ) for problem sizes from 1116 to 1216 500 triples (α, β, ρ) for problem sizes from 1116 to 1516 800 triples (α, β, ρ) for problem sizes from 1116 to 1816 1000 triples (α, β, ρ) for problem sizes from 1116 to 2016

8 7

Run time, s

6 5

4 3 2 1 0 1116

1216

1316

1416

1516

1616

1716

1816

1916

Number of CES variants (problem size)

2016

Figure 6.7: Run time of the program for problems of various sizes devices have been synchronized (line 23), the solution found is copied to the host from the device on which it was obtained (line 26). If no device has found a solution within the allotted number of iterations, the attempt is repeated (line 24), as in the single-GPU version (see Listing 5.3, lines 16–20). In terms of the classification given in Section 3.1, this scheme is static: from the very beginning, each device receives its own fixed set of independent colonies, and no work is redistributed between devices. For the problem under consideration, this causes no losses due to load imbalance, since the search terminates not after all work has been done but at the first success of any colony.

6.5.2

Experimental Platforms

The specifications of the GPUs used are given in Table 6.3. The peak single-precision performance was computed as twice the product of the number of CUDA cores and the clock rate. Table 6.3: Specifications of the GPUs used in the experiments Specification

Tesla K20m

TITAN Xp

Tesla V100-SXM2

Microarchitecture Compute Capability Number of streaming multiprocessors CUDA cores per multiprocessor Total CUDA cores Clock rate, MHz Global memory, GB Peak performance (FP32), TFLOPS

Kepler 3.5 13 192 2496 705 5 3.5

Pascal 6.1 30 128 3840 1582 12 12.1

Volta 7.0 80 64 5120 1530 16 15.7

Number of GPUs in the compute node Total number of colonies using all GPUs

3 7488

8 30720

4 20480

79

7

Multi-constant: 200 triples (α, β, ρ) for problem sizes from 1116 to 1216 Multi-constant: 1000 triples (α, β, ρ) for problem sizes from 1116 to 2016 Random: α ∈ (0.1, 0.3], β ∈ (0.3, 0.7], ρ ∈ (0.8, 1.0] Constant: α = 0.2, β = 0.5, ρ = 0.9

6 5

Run time, s

4 3 2 1 0 1116

1216

1316

1416

1516

1616

1716

1816

Number of CES variants (problem size)

1916

2016

Figure 6.8: Run time of the program for problems of various sizes: comparison of the three tuning approaches

6.5.3

Comparison of GPU Architectures

Table 6.4 gives the time to find a feasible solution on a single accelerator, averaged over all problem sizes (from 1116 to 2016 ), for the different approaches to parameter tuning. For the TITAN Xp, results are available for only two parameter settings. The search time on the Tesla V100 was 5.8–9.1 times shorter than on the Tesla K20m and 2.5–2.9 times shorter than on the TITAN Xp. These ratios exceed the ratios of the peak performance of the accelerators (4.5 and 1.3, respectively; see Table 6.3). The reason is that a more powerful accelerator not only performs the computations of each colony faster but also allows more colonies to run simultaneously: 5120 on the V100 versus 3840 on the TITAN Xp and 2496 on the K20m. As the number of colonies grows, the probability that at least one of them finds a feasible solution before the allotted number of iterations is exhausted increases, and the number of repeated attempts decreases: for example, for the intervals α ∈ (0, 2], β ∈ (0, 2], ρ ∈ (0, 1], the average number of attempts was 3.7 on the K20m, 2.7 on the TITAN Xp, and 2.3 on the V100. At the same time, the relative efficiency of the tuning approaches does not depend on the architecture. On the K20m, as on the V100, narrowing the parameter intervals in the random approach successively reduces the search time (from 339.5 to 21.3 s, i.e., by a factor of 16), and the best result is achieved by the constant approach with the values α = 0.2, β = 0.5, ρ = 0.9 found in Section 6.4: compared with the original values α = 0.4, β = 0.6, ρ = 0.9, the search time on the K20m decreased by a factor of 34 (from 347.5 to 10.1 s). Thus, parameter values obtained by tuning on one architecture remain effective on another.

6.5.4

Scaling the Parallel Program to Multiple GPUs

To assess scalability, we used the random approach with the original intervals α ∈ (0, 2], β ∈ (0, 2], ρ ∈ (0, 1], for which results are available on all three platforms, as well as the constant approach with α = 0.4, β = 0.6, ρ = 0.9 on the Tesla K20m platform. The results are shown in Figs. 6.9 and 6.10 and in Table 6.5, where T1 and TN are the average search times on one and on N 80

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

F i n d F i r s t S o l u t i o n M u l t i G P U () { for ( d = 1; d <= numDevices ; d ++) { /* prepare each device */ cudaSetDevice ( d ) ; cud aMemcpyH toD ( inData , oprData ) ; initRandomStates < < < blocksPerGrid , threadsPerBlock > > >( seed * d ) ; SET < < < blocksPerGrid , threadsPerBlock > > >( TUN_APP ) ; } do { for ( d = 1; d <= numDevices ; d ++) { /* a s y n c h r o n o u s ACO launch */ cudaSetDevice ( d ) ; ACO < < < blocksPerGrid , threadsPerBlock , 0 , streamKernel [ d ] > > >() ; cu da Ev e nt Re co r d ( kernelDone [ d ] , streamKernel [ d ]) ; } finished = 0; bestDevice = -1; while ( finished < numDevices ) { /* round - robin polling of the devices */ for ( d = 1; d <= numDevices ; d ++) { if (! isFinished [ d ] && cud aEventQu ery ( kernelDone [ d ]) == cudaSuccess ) { isFinished [ d ] = true ; finished ++; cud aMemcpyD toH ( isFound , d ) ; if ( bestDevice == -1 && isFound != 0) { bestDevice = d ; /* stop the search on all devices */ for ( k = 1; k <= numDevices ; k ++) c u d a M e m c p y H t o D A s y n c ( isFound , k , streamStop [ k ]) ; } } } } c u d a D e v i c e S y n c h r o n i z e A l l () ; } while ( bestDevice == -1) ; cudaSetDevice ( bestDevice ) ; cuda MemcpyDt oH ( Wfirst , threadIdx ) ; }

Listing 6.2: Pseudocode of the search for a feasible solution on multiple GPUs (host)

accelerators, S = T1 /TN is the speedup, and E = S/N is the efficiency. Figures 6.9 and 6.10 refer to the random approach with the original parameter intervals; because of the large spread of values, the search time in Fig. 6.9 is shown on a logarithmic scale. In all cases, the speedup was at least linear and, in most cases, superlinear: the efficiency ranged from 1.03 to 1.51. The largest speedup was obtained on the node with eight TITAN Xp accelerators: the average search time decreased from 125.9 to 10.6 s, i.e., by a factor of 11.9. For individual problem sizes, the speedup was even higher: for example, for the problem size 1116 , it was 28 on eight TITAN Xp and 11.6 on four V100 (Fig. 6.10). The superlinear speedup is not a measurement artifact but follows from the nature of the problem. As noted in Section 3.1.5, if the search stops when the first solution is found, the parallel algorithm may perform less work than the sequential one, which corresponds to a search overhead factor R < 1 (expression (3.1)). In this case, the run time is determined not by the total amount of computation but by the time to the first success among all simultaneously running colonies, i.e., by the minimum of random variables whose number is proportional to the number of accelerators. If the time a single colony needs to find a solution varies strongly from run to run (such distributions are called heavy-tailed), then, alongside long runs, there is a noticeable fraction of very short ones, and increasing the number of independent colonies reduces the expected search time faster than linearly [190]. A similar effect is well known for parallel evolutionary and other population-based algorithms [191]. In practice, it manifests itself as a reduction in the number of repeated attempts: when all accelerators of the node are used, a feasible solution is usually found at the very first kernel launch (Table 6.5), whereas on a single GPU an average of 2.3 to 4.3 attempts were needed in the cases considered. Note also that the node with eight TITAN Xp accelerators (30 720 colonies) and the node with four V100 accelerators (20 480 colonies) achieved practically the same search time (10.6 and 10.0 s), even though on a single GPU the V100 is 2.5 times faster than the TITAN Xp. Hence, for the problem under consideration, the lower performance of an individual accelerator can be largely compensated for by increasing the total number of simultaneously running colonies. Combining multiple GPUs with parameter tuning gives an additional gain. On the node with eight TITAN Xp, using the narrowed intervals of the random approach, the average search time 81

Table 6.4: Average time to find a feasible solution on a single GPU, s Approach and parameter values

K20m

TITAN Xp

V100

K20m / V100

Random: α ∈ (0, 2], β ∈ (0, 2], ρ ∈ (0, 1]

339.5

125.9

51.1

6.6

Random: α ∈ (0, 1.2], β ∈ (0.1, 1.2], ρ ∈ (0.2, 1]

108.9

–

13.9

7.9

Random: α ∈ (0, 0.8], β ∈ (0.1, 1.1], ρ ∈ (0.3, 1]

57.6

–

7.7

7.5

Random: α ∈ (0, 0.6], β ∈ (0.2, 1.1], ρ ∈ (0.4, 1]

44.7

–

5.5

8.1

Random: α ∈ (0, 0.5], β ∈ (0.2, 0.9], ρ ∈ (0.5, 1]

28.7

–

3.4

8.4

Random: α ∈ (0, 0.4], β ∈ (0.2, 0.8], ρ ∈ (0.6, 1]

21.3

–

2.7

7.7

Constant: α = 0.4, β = 0.6, ρ = 0.9

347.5

111.2

38.1

9.1

Constant: α = 0.2, β = 0.5, ρ = 0.9

10.1

–

1.5

6.5

Multi-constant: 500 triples for problem sizes from 1116 to 1516

18.3

–

3.1

5.8

Table 6.5: Speedup and efficiency with several GPUs GPU, parameters

N

T1 , s

TN , s

S

E

Attempts (1 / N GPUs)

K20m, random K20m, constant TITAN Xp, random V100, random

3 3 8 4

339.5 347.5 125.9 51.1

110.3 76.7 10.6 10.0

3.1 4.5 11.9 5.1

1.03 1.51 1.48 1.28

3.7 / 1.7 4.3 / 1.7 2.7 / 1.0 2.3 / 1.1

ranged from 1.3 to 2.7 s: for example, 1.3 s for the intervals α ∈ (0, 0.5], β ∈ (0.2, 0.9], ρ ∈ (0.5, 1] versus 3.4 s on a single V100, and 2.7 s for the intervals α ∈ (0, 1.2], β ∈ (0.1, 1.2], ρ ∈ (0.2, 1] versus 13.9 s. With the multi-constant approach using 500 parameter triples, the search time on eight TITAN Xp was 1.9 s versus 3.1 s on a single V100. The study on GPUs of three generations showed that the parameter values found remain effective regardless of the GPU architecture, and that using several accelerators within a single compute node provides at least linear speedup of the search for a feasible solution (up to 11.9 times on eight TITAN Xp accelerators) without changing the compute kernel. These results show that the proposed parallel algorithm for finding a feasible solution scales well to multiple GPUs without any changes to the compute kernel. Since the interaction between devices reduces to sending a flag when a solution is found and copying the solution to the host, the communication overhead is negligible, and increasing the number of accelerators in a node is a simple and effective way of reducing the search time for large problems.

82

Tesla K20m, 1 GPU Tesla K20m, 3 GPUs

10000

TITAN Xp, 1 GPU TITAN Xp, 8 GPUs

Tesla V100, 1 GPU Tesla V100, 4 GPUs

Run time, s

1000

100

10

1 1116

1216

1316

1416

1516

1616

1716

1816

1916

Number of CES variants (problem size)

Figure 6.9: Time to find a feasible solution on one and on several GPUs

35

Tesla K20m, 3 GPUs TITAN Xp, 8 GPUs Tesla V100, 4 GPUs

30 25

Speedup

20 15

10 5 0

1116

1216

1316

1416

1516

1616

1716

1816

Number of CES variants (problem size)

1916

2016

Figure 6.10: Speedup of the search for a feasible solution on several GPUs

83

2016

Conclusion The problem of optimal equipment selection for chemical-engineering systems of multiproduct batch plants has been formulated (Chapter 1). A sequential algorithm based on the branch-and-bound method has been developed for the optimal CES equipment selection problem (Chapter 2). The algorithm has been implemented in C++. The performance of individual elements of the algorithm has been studied, and on this basis the algorithm has been optimized to increase its speed (Chapter 2). The algorithm has been parallelized. For shared-memory systems, a program has been developed in C++ using the OpenMP standard (Section 3.2). The speedup and processor efficiency of this program have been studied. The study showed that the parallel program achieves near-linear speedup and fairly high scalability, with a core efficiency of 0.9–0.95. For distributed-memory systems, a program has been developed in C++ using the MPI library, implementing the algorithm according to the master–worker paradigm (Section 3.3). The study of speedup and efficiency showed that the parallel program achieves near-linear speedup and fairly high scalability. However, the master process becomes a bottleneck of the program when it serves a large number of workers, which reduces the speedup and efficiency of the program. For heterogeneous computing systems equipped with GPUs, a parallel program has been developed in C using CUDA, taking into account the technical features of GPUs (Chapter 4). Two ways of traversing tree structures have been implemented: recursive and iterative. Two ways of optimizing the algorithm have also been proposed: efficient use of high-performance shared memory and reduction of branch divergence. The efficiency of using GPUs for the problem under consideration has been evaluated. Although the program accelerates the solution of the problem considered in this book, a direct implementation of branch-and-bound on the GPU does not yield a large speedup (by hundreds or thousands of times). A hybrid metaheuristic algorithm (ACO+SA) combining two metaheuristics, Ant Colony Optimization (ACO) and Simulated Annealing (SA), has been developed (Chapter 5). The hybrid algorithm works in two stages: first, a feasible initial solution (a fully operable CES variant that produces the given product range within the given time) is found by ACO, and then this solution is improved by SA. In developing the algorithm, an original way of computing the fitness function for ant colony optimization, on the basis of which the pheromone concentration is updated, has been proposed, together with a way of defining the heuristic information. For simulated annealing, an original way of searching for a candidate solution in the neighborhood of the current feasible solution has been proposed, as well as a way of determining the maximum and minimum temperatures of the cooling schedule. The algorithm has been implemented in C using the CUDA programming model. The program includes code executed on the CPU, which implements only the sequential parts of the algorithm, prepares data and copies them to the GPU, sets the kernel parameters, and launches the kernels. The efficiency of using GPUs for solving the problem with the hybrid algorithm has been evaluated. The parallel implementation of the hybrid metaheuristic algorithm (ACO+SA) yields high-quality solutions for the test problem that coincide with the globally optimal values obtained by the deterministic branch-and-bound algorithm, but much faster: in a few tens of seconds instead of a few tens of hours. However, when the problem size exceeds 84

1116 , the performance drops considerably. This is because, as the problem size grows, finding a feasible initial solution with ACO takes more and more time. Problems larger than 1216 were not solved, which motivated the transition to tuning the parameters of the algorithm. Before ACO is started, a number of parameters that affect its operation must be set (tuned) (Chapter 6). A family of methods for tuning the parameters α, β, and ρ of the algorithm has been developed: random, constant, and multi-constant. Owing to modern heterogeneous CPU–GPU computing systems, a large number of computational experiments could be performed in relatively short time, and their results were then processed statistically. By analyzing how often solutions fall into the ranges of the algorithm parameters, it was possible to set the bounds of these ranges to values that yield the largest number of solutions and to accelerate the algorithm considerably. Several successive series of experiments with the random approach yielded values of the constants for the constant approach that provide the highest performance. To initialize the algorithm parameters in the multi-constant approach, the values obtained in the first series of experiments with the random approach were used. The study of the multi-constant approach showed that parameters tuned for small problems can be used to solve large problems. Although this approach is not the best one in terms of algorithm speed, it considerably reduces the time spent on searching for the best algorithm parameters. The parameter values found for the metaheuristic algorithm reduced the average time needed to find a feasible initial solution for problems of size 1116 to 2016 by more than an order of magnitude: from ≈ 38.1 s for α = 0.4, β = 0.6, ρ = 0.9 to ≈ 1.5 s for α = 0.2, β = 0.5, ρ = 0.9. The multi-constant approach provides a search time of ≈ 2.9–3.2 s while reducing the time spent on tuning the parameters by a factor of ≈ 26. The program has additionally been studied on GPUs of three generations (Tesla K20m, TITAN Xp, and Tesla V100) and on compute nodes with several GPUs (Section 6.5). The parameter values found by tuning on one architecture remain effective on the others: on the Tesla K20m, they reduce the average time needed to find a feasible solution by a factor of 34 (from 347.5 to 10.1 s). Since the ant colonies are independent, the algorithm extends to several GPUs without any changes to the kernel function. The speedup obtained is at least linear and, in most cases, superlinear: up to 11.9 on eight TITAN Xp accelerators (from 125.9 to 10.6 s). The reason is that the search time is determined by the first success among all simultaneously running colonies, and the number of repeated attempts decreases as the number of colonies grows.

85

Bibliography [1] N. Fujimoto and Sh. Tsutsui. A highly-parallel TSP solver for a GPU computing platform. In Numerical Methods and Applications, pages 264–271. Springer, 2011. doi: 10.1007/978-3 -642-18466-6_31. [2] Michał Czapiński. An effective parallel multistart tabu search for quadratic assignment problem on CUDA platform. Journal of Parallel and Distributed Computing, 73(11): 1461–1468, 2013. doi: 10.1016/j.jpdc.2012.07.014. [3] Vincent Boyer, Didier El Baz, and Moussa Elkihel. Solving knapsack problems on GPU. Computers & Operations Research, 39(1):42–47, 2012. doi: 10.1016/j.cor.2011.03.014. [4] Robert Pellerin, Nathalie Perrier, and François Berthaut. A survey of hybrid metaheuristics for the resource-constrained project scheduling problem. European Journal of Operational Research, 280(2):395–416, 2020. doi: 10.1016/j.ejor.2019.01.063. [5] M. Gavanelli and T. Mancini. Experimental Evaluation of Algorithms for Solving Problems with Combinatorial Explosion: RCRA 2007. Elsevier Incorporated, 2009. Special Issue. [6] Gerhard J. Woeginger. Exact Algorithms for NP-Hard Problems: A Survey, pages 185– 207. Springer Berlin Heidelberg, Berlin, Heidelberg, 2003. ISBN 978-3-540-36478-8. doi: 10.1007/3-540-36478-1_17. [7] Janez Žerovnik. Heuristics for NP-hard optimization problems – simpler is better!? Logistics & Sustainable Transport, 6(1):1–10, 2015. doi: 10.1515/jlst-2015-0006. [8] E.-G. Talbi. Metaheuristics: from design to implementation. John Wiley & Sons, 2009. ISBN 978-0-470-27858-1. [9] Jürg Nievergelt. Exhaustive Search, Combinatorial Optimization and Enumeration: Exploring the Potential of Raw Computing Power. In Václav Hlaváč, Keith G. Jeffery, and Jiří Wiedermann, editors, SOFSEM 2000: Theory and Practice of Informatics, pages 18–35, Berlin, Heidelberg, 2000. Springer Berlin Heidelberg. ISBN 978-3-540-44411-4. doi: 10.1007/3-540-44411-4_2. [10] Caio Paziani Tomazella and Marcelo Seido Nagano. A comprehensive review of Branchand-Bound algorithms: Guidelines and directions for further research on the flowshop scheduling problem. Expert Systems with Applications, 158:113556, 2020. doi: 10.1016/j.es wa.2020.113556. [11] Ketan Date and Rakesh Nagi. GPU-accelerated Hungarian algorithms for the Linear Assignment Problem. Parallel Computing, 57:52–72, 2016. doi: 10.1016/j.parco.2016.05.012. [12] Paul Bouman, Niels Agatz, and Marie Schmidt. Dynamic programming approaches for the traveling salesman problem with drone. Networks, 72(4):528–542, 2018. doi: 10.1002/net.21864.

86

[13] Rafael Martí and Gerhard Reinelt. Heuristic Methods, pages 17–40. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011. ISBN 978-3-642-16729-4. doi: 10.1007/978-3-642-167 29-4_2. [14] A. V. Panteleev. Metaheuristic Algorithms for Global Extremum Search. MAI-Print, Moscow, 2009. ISBN 978-5-7035-2115-1. In Russian. [15] Xin-She Yang. Nature-inspired metaheuristic algorithms. Luniver press, 2010. [16] N. M. Ershov and N. N. Popova. Natural Models of Parallel Computations. Computer Research and Modeling, 7(3):781–785, 2015. doi: 10.20537/2076-7633-2015-7-3-781-785. In Russian. [17] Fred Glover and Manuel Laguna. Tabu Search. In Handbook of Combinatorial Optimization, pages 3261–3362. Springer New York, New York, NY, 2013. doi: 10.1007/978-1-4419-799 7-1_17. [18] Emile Aarts, Jan Korst, and Wil Michiels. Simulated Annealing. In Search Methodologies, pages 265–285. Springer US, Boston, MA, 2014. doi: 10.1007/978-1-4614-6940-7_10. [19] Marco Dorigo and Christian Blum. Ant Colony Optimization Theory: A Survey. Theoretical Computer Science, 344(2–3):243–278, 2005. doi: 10.1016/j.tcs.2005.05.020. [20] Kim-Fung Man, Kit Sang Tang, and Sam Kwong. Genetic Algorithms. Advanced Textbooks in Control and Signal Processing. Springer London, London, 1999. ISBN 978-1-85233-072-9. doi: 10.1007/978-1-4471-0577-0. [21] Maurice Clerc. Particle Swarm Optimization. ISTE / John Wiley & Sons, 2013. ISBN 978-1-118-61397-9. [22] Dervis Karaboga, Beyza Gorkemli, Celal Ozturk, and Nurhan Karaboga. A comprehensive survey: artificial bee colony (ABC) algorithm and applications. Artificial Intelligence Review, 42(1):21–57, 2014. doi: 10.1007/s10462-012-9328-0. [23] Davide Anghinolfi and Massimo Paolucci. Simulated Annealing as an Intensification Component in Hybrid Population-based Metaheuristics. In Simulated Annealing, pages 1–28. InTech, 2008. doi: 10.5772/5570. [24] Christian Blum, Jakob Puchinger, Günther R. Raidl, and Andrea Roli. Hybrid metaheuristics in combinatorial optimization: A survey. Applied Soft Computing, 11(6):4135–4151, 2011. doi: 10.1016/j.asoc.2011.02.032. [25] Mauro Birattari. Tuning Metaheuristics, volume 197 of Studies in Computational Intelligence. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009. ISBN 978-3-642-00482-7. doi: 10.1007/978-3-642-00483-4. [26] Susheel Kumar Joshi and Jagdish Chand Bansal. Parameter Tuning for Meta-Heuristics. Knowledge-Based Systems, 189:105094, 2020. doi: 10.1016/j.knosys.2019.105094. [27] Changwu Huang, Yuanxiang Li, and Xin Yao. A Survey of Automatic Parameter Tuning Methods for Metaheuristics. IEEE Transactions on Evolutionary Computation, 24(2): 201–216, 2020. doi: 10.1109/TEVC.2019.2921598. [28] V. V. Voevodin and Vl. V. Voevodin. Parallel Computing. BHV-Peterburg, Saint Petersburg, 2002. ISBN 5-94157-160-7. In Russian.

87

[29] Ananth Grama, Anshul Gupta, George Karypis, and Vipin Kumar. Introduction to Parallel Computing, Design and Analysis of Algorithms. Pearson Education, second edition, 2003. ISBN 9788131708071. [30] Fayez Gebali. Algorithms and parallel computing, volume 84. John Wiley & Sons, 2011. [31] R. G. Strongin, V. P. Gergel, V. A. Grishagin, and K. A. Barkalov. Parallel Computations in Global Optimization Problems. Moscow University Press, Moscow, 2013. ISBN 978-5211-06479-9. In Russian. [32] David Padua, editor. Encyclopedia of Parallel Computing. Springer US, Boston, MA, 2011. ISBN 978-0-387-09765-7. doi: 10.1007/978-0-387-09766-4. [33] Imen Chakroun, Nordine Melab, Mohand Mezmaz, and Daniel Tuyttens. Combining multi-core and GPU computing for solving combinatorial optimization problems. Journal of Parallel and Distributed Computing, 73(12):1563–1577, 2013. doi: 10.1016/j.jpdc.2013. 07.023. [34] A. A. Kolpakov and Yu. A. Kropotov. Improving the Performance of Heterogeneous Computer Data Processing Systems. Direkt-Media, Moscow; Berlin, 2019. ISBN 978-54475-8156-5. In Russian. [35] E. V. Semenov, A. V. Starchenko, E. A. Danilkin, and S. A. Prokhanov. Choosing a Parallelization Technology for the Numerical Solution of the Convection–Diffusion Equation on a Hybrid Multiprocessor Computing System. Tomsk State University Journal of Control and Computer Science, (51):111–120, 2020. doi: 10.17223/19988605/51/13. In Russian. [36] Davor Sluga, Tomaz Curk, Blaz Zupan, and Uros Lotric. Heterogeneous computing architecture for fast detection of SNP-SNP interactions. BMC bioinformatics, 15(1):216, 2014. [37] Gabriela Sandoval, Daniel Espinoza, Nicolas Figueroa, and Juan A. Asenjo. Optimization of a biotechnological multiproduct batch plant design for the manufacture of four different products: A real case scenario. Biotechnology and Bioengineering, 114(6):1252–1263, 2017. doi: 10.1002/bit.26260. [38] A. B. Borisenko and S. V. Karpushkin. Hierarchy of processing equipment configuration design problems for multiproduct chemical plants. Journal of Computer and Systems Sciences International, 53(3):410–419, 2014. doi: 10.1134/S1064230714030046. Translated from Izvestiya Rossiiskoi Akademii Nauk. Teoriya i Sistemy Upravleniya, 2014, no. 3, pp. 113–123. [39] S. V. Karpushkin, M. N. Krasnyanskiy, and A. B. Borisenko. Retrofitting the ChemicalEngineering Systems of an Operating Multiproduct Plant When Production Plans Change. Khimicheskaya Promyshlennost’ Segodnya [Chemical Industry Today], (1):26–32, 2012. In Russian. [40] V. V. Kafarov and V. V. Makarov. Flexible Automated Production Systems in the Chemical Industry. Khimiya, Moscow, 1990. ISBN 5-7245-0236-4. In Russian. [41] R. Fletcher, J. A. J. Hall, and W. R. Johns. Flexible retrofit design of multiproduct batch plants. Computers & Chemical Engineering, 15(12):843–852, 1991. doi: 10.1016/0098-135 4(91)80029-u. [42] E. N. Malygin, S. V. Karpushkin, and A. B. Borisenko. Design of Multiproduct Chemical Plants: Determining the Equipment Configuration of Chemical-Engineering Flowsheets. 88

Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 8(2):272–283, 2002. In Russian. [43] E. N. Malygin, A. B. Borisenko, and S. V. Karpushkin. A Method for Determining the Equipment Configuration of Multiproduct Chemical-Engineering Systems. Khimicheskaya Promyshlennost’ Segodnya [Chemical Industry Today], (5):43–50, 2003. In Russian. [44] S. V. Karpushkin, M. N. Krasnyanskiy, and A. B. Borisenko. A Method for Assessing the Efficiency of Equipment Configurations of Chemical-Engineering Systems of a Multiproduct Plant. Informatsionnye Sistemy i Tekhnologii [Information Systems and Technologies], (5(67)):96–105, 2011. In Russian. [45] E. N. Malygin, S. V. Karpushkin, and P. G. Mikhailova. Selection of Auxiliary Process Equipment for Multiproduct Chemical Plants. Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 17(2): 483–487, 2011. In Russian. [46] A. B. Borisenko, A. V. Antonenko, A. V. Osovskiy, and O. A. Filimonova. A System for the Computer-Aided Selection of Auxiliary Equipment for Multiproduct Chemical Plants. Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 18(3):569–572, 2012. In Russian. [47] Youness El Hamzaoui, J. A. Hernandez, M. A. Cruz-Chavez, and A. Bassam. Search for Optimal Design of Multiproduct Batch Plants under Uncertain Demand using Gaussian Process Modeling Solved by Heuristics Methods. Chemical Product and Process Modeling, 5(1), 2010. doi: 10.2202/1934-2659.1426. Article 8. [48] Antonin Ponsich, Catherine Azzaro-Pantel, Serge Domenech, and Luc Pibouleau. Mixedinteger nonlinear programming optimization strategies for batch plant design problems. Industrial and Engineering Chemistry Research, 46(3):854–863, 2007. doi: 10.1021/ie0607 33d. [49] Yunfei Chu and Fengqi You. Integrated Planning, Scheduling, and Dynamic Optimization for Batch Processes: MINLP Model Formulation and Efficient Solution Methods via Surrogate Modeling. Industrial & Engineering Chemistry Research, 53(34):13391–13411, 2014. doi: 10.1021/ie501986d. [50] G. Sandoval, D. Espinoza, N. Figueroa, and J. A. Asenjo. MILP reformulations for the design of biotechnological multi-product batch plants using continuous equipment sizes and discrete host selection. Computers & Chemical Engineering, 84:1–11, 2016. doi: 10.1016/j.compchemeng.2015.08.001. [51] Yanina Fumero, Gabriela Corsano, and Jorge M. Montagna. A Mixed Integer Linear Programming model for simultaneous design and scheduling of flowshop plants. Applied Mathematical Modelling, 37(4):1652–1664, 2013. doi: 10.1016/j.apm.2012.04.043. [52] Xiaoxia Lin and C. A. Floudas. Design, synthesis and scheduling of multipurpose batch plants via an effective continuous-time formulation. In Computers and Chemical Engineering, volume 25, pages 665–674. Elsevier, 2001. doi: 10.1016/S0098-1354(01)00663-9. [53] Esmael R. Seid and Thokozani Majozi. Design and Synthesis of Multipurpose Batch Plants Using a Robust Scheduling Platform. Industrial & Engineering Chemistry Research, 52 (46):16301–16313, 2013. doi: 10.1021/ie4022495. [54] Andrey Borisenko, Philipp Kegel, and Sergei Gorlatch. Optimal Design of Multi-product Batch Plants Using a Parallel Branch-and-Bound Method. In Victor Malyshkin, editor, 89

Parallel Computing Technologies (PaCT 2011), volume 6873 of Lecture Notes in Computer Science, pages 417–430, Berlin, Heidelberg, 2011. Springer. doi: 10.1007/978-3-642-23178 -0_36. [55] A. B. Borisenko and S. V. Karpushkin. Application of the Branch-and-Bound Method to the Optimal Equipment Selection for Chemical-Engineering Systems. Vychislitel’nye Tekhnologii [Computational Technologies], 17(1):35–43, 2012. In Russian. [56] João Luís de Miranda. The design and scheduling of chemical batch processes: Computational complexity studies. Computers & Chemical Engineering, 121:367–374, 2019. doi: 10.1016/j.compchemeng.2018.11.011. [57] Andrey Borisenko and Sergei Gorlatch. Optimal Batch Plants Design on Parallel Systems: A Comparative Study. In 2019 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 549–558, 2019. ISBN 978-1-7281-3510-6. doi: 10.1109/IPDPSW.2019.00097. [58] E. N. Malygin, S. V. Karpushkin, and A. B. Borisenko. A mathematical model of the functioning of multiproduct chemical engineering systems. Theoretical Foundations of Chemical Engineering, 39(4):429–439, 2005. doi: 10.1007/s11236-005-0099-z. Translated from Teoreticheskie Osnovy Khimicheskoi Tekhnologii, 2005, vol. 39, no. 4, pp. 455–465. [59] D. W. T. Rippin. Design and Operation of Multiproduct and Multipurpose Batch Chemical Plants. – Analysis of Problem Structure. Computers & Chemical Engineering, 7(4):463–481, 1983. doi: 10.1016/0098-1354(83)80023-4. [60] A. P. Barbosa-Póvoa. A Critical Review on the Design and Retrofit of Batch Plants. Computers & Chemical Engineering, 31(7):833–855, 2007. doi: 10.1016/j.compchemeng.20 06.08.003. [61] G. R. Kocis and I. E. Grossmann. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Industrial & Engineering Chemistry Research, 27(8):1407–1421, 1988. doi: 10.1021/ie00080a013. [62] D. E. Ravemark and D. W. T. Rippin. Predesigning a Multiproduct Batch Plant by Mathematical Programming. In G. V. Reklaitis, A. K. Sunol, D. W. T. Rippin, and Ö. Hortaçsu, editors, Batch Processing Systems Engineering, volume 143 of NATO ASI Series, pages 114–149. Springer, Berlin, Heidelberg, 1996. doi: 10.1007/978-3-642-60972-5 _7. [63] J. M. Montagna, A. R. Vecchietti, O. A. Iribarren, J. M. Pinto, and J. A. Asenjo. Optimal Design of Protein Production Plants with Time and Size Factor Process Models. Biotechnology Progress, 16(2):228–237, 2000. doi: 10.1021/bp990145d. [64] M. S. Moreno, O. A. Iribarren, and J. M. Montagna. Optimal Design of Multiproduct Batch Plants Considering Duplication of Units in Series. Chemical Engineering Research and Design, 87:1497–1508, 2009. doi: 10.1016/j.cherd.2009.04.014. [65] E. N. Malygin and S. V. Karpushkin. Design of Multiproduct Chemical Plants: Determining the Cycle Times for Processing Product Batches. Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 5(2): 201–212, 1999. In Russian. [66] S. V. Karpushkin. Equipment Selection for Multiproduct Chemical Plants. Mashinostroenie1, Moscow, 2006. ISBN 5-94275-262-1. In Russian.

90

[67] S. V. Karpushkin. A Methodology for the Equipment Design of Multiproduct Chemical Plants. Dr. sci. (eng.) dissertation, Tambov State Technical University, Tambov, 2007. In Russian. [68] Gilles Brassard and Paul Bratley. Fundamentals of Algorithmics. Prentice Hall, 1996. [69] David R. Morrison, Sheldon H. Jacobson, Jason J. Sauppe, and Edward C. Sewell. Branchand-bound algorithms: A survey of recent advances in searching, branching, and pruning. Discrete Optimization, 19:79–102, 2016. doi: 10.1016/j.disopt.2016.01.005. [70] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algorithms. MIT Press, Cambridge, MA, 3rd edition, 2009. ISBN 978-0-262-03384-8. [71] D. E. Knuth. The Art of Computer Programming, Vol. 3: Sorting and Searching. AddisonWesley, Reading, MA, 2nd edition, 1998. ISBN 978-0-201-89685-5. [72] Sven Leyffer, Jeff Linderoth, James Luedtke, Andrew Miller, and Todd Munson. Applications and algorithms for mixed integer nonlinear programming. Journal of Physics: Conference Series, 180:012014, 2009. doi: 10.1088/1742-6596/180/1/012014. [73] Pierre Bonami, Mustafa Kilinç, and Jeff Linderoth. Algorithms and software for convex mixed integer nonlinear programs. In Mixed Integer Nonlinear Programming, pages 1–39. Springer, 2012. doi: 10.1007/978-1-4614-1927-3_1. [74] Samuel Burer and Adam N. Letchford. Non-convex mixed-integer nonlinear programming: A survey. Surveys in Operations Research and Management Science, 17(2):97–106, 2012. doi: 10.1016/j.sorms.2012.08.001. [75] M. A. Posypkin and I. Kh. Sigal. Investigation of Algorithms for Parallel Computations in Knapsack-Type Discrete Optimization Problems. Computational Mathematics and Mathematical Physics, 45(10):1735–1742, 2005. Translated from Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki, 2005, vol. 45, no. 10, pp. 1801–1809. [76] A. B. Borisenko, D. V. Kutuzov, and A. V. Osovskiy. Application of Parallel Computing to the Calculation of the Equipment Configuration of Chemical-Engineering Systems. Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 17(2):493–496, 2011. In Russian. [77] N. A. Ezhova and L. B. Sokolinsky. Parallel Computational Model for Multiprocessor Systems with Distributed Memory. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering, 7(2):32–49, 2018. doi: 10.14529/c mse180203. In Russian. [78] Ananth Grama and Vipin Kumar. State of the Art in Parallel Search Techniques for Discrete Optimization Problems. IEEE Transactions on Knowledge and Data Engineering, 11(1):28–35, 1999. doi: 10.1109/69.755612. [79] Bernard Gendron and Teodor Gabriel Crainic. Parallel Branch-and-Bound Algorithms: Survey and Synthesis. Operations Research, 42(6):1042–1066, 1994. doi: 10.1287/opre.42.6. 1042. [80] Vipin Kumar, Ananth Y. Grama, and Nageshwara Rao Vempaty. Scalable Load Balancing Techniques for Parallel Computers. Journal of Parallel and Distributed Computing, 22(1): 60–79, 1994. [81] N. E. Timoshevskaya. Parallel Tree Traversal Methods. Matematicheskoe Modelirovanie [Mathematical Modeling], 16(1):105–114, 2004. In Russian. 91

[82] A. B. Borisenko, S. V. Karpushkin, and A. O. Glebov. A Parallel Algorithm for Solving the Three-Dimensional Unsteady Heat Equation with an Explicit Difference Scheme. Vestnik Tambovskogo Gosudarstvennogo Tekhnicheskogo Universiteta [Transactions of Tambov State Technical University], 16(3):573–577, 2010. In Russian. [83] Teodor Gabriel Crainic, Bertrand Le Cun, and Catherine Roucairol. Parallel Branch-andBound Algorithms. In El-Ghazali Talbi, editor, Parallel Combinatorial Optimization, pages 1–28. Wiley, Hoboken, 2006. doi: 10.1002/9780470053928.ch1. [84] Wei Shu and L. V. Kale. A Dynamic Scheduling Strategy for the Chare-Kernel System. In Proceedings of the 1989 ACM/IEEE Conference on Supercomputing, pages 389–398, 1989. [85] V. Nageshwara Rao and Vipin Kumar. Parallel Depth First Search. Part I. Implementation. International Journal of Parallel Programming, 16(6):479–499, 1987. [86] Alexander Reinefeld and Volker Schnecke. Work-Load Balancing in Highly Parallel DepthFirst Search. In Proceedings of the Scalable High-Performance Computing Conference, pages 773–780, 1994. [87] Raphael Finkel and Udi Manber. DIB – A Distributed Implementation of Backtracking. ACM Transactions on Programming Languages and Systems, 9(2):235–256, 1987. [88] Peter Sanders. Better Algorithms for Parallel Backtracking. In Parallel Algorithms for Irregularly Structured Problems, volume 980 of Lecture Notes in Computer Science, pages 333–347. Springer, 1995. [89] Robert D. Blumofe and Charles E. Leiserson. Scheduling Multithreaded Computations by Work Stealing. Journal of the ACM, 46(5):720–748, 1999. [90] V. Nageshwara Rao and Vipin Kumar. On the Efficiency of Parallel Backtracking. IEEE Transactions on Parallel and Distributed Systems, 4(4):427–437, 1993. [91] Ten-Hwang Lai and Sartaj Sahni. Anomalies in Parallel Branch-and-Bound Algorithms. Communications of the ACM, 27(6):594–602, 1984. [92] N. E. Timoshevskaya. Development and Study of Parallel Combinatorial Algorithms. Prikladnaya Diskretnaya Matematika [Applied Discrete Mathematics], (2(4)):96–103, 2009. In Russian. [93] Mohamed Esseghir Lalami and Didier El-Baz. GPU Implementation of the Branch and Bound Method for Knapsack Problems. In 2012 IEEE 26th International Parallel and Distributed Processing Symposium Workshops & PhD Forum, pages 1763–1771, 2012. doi: 10.1109/IPDPSW.2012.219. [94] Andrey Borisenko, Michael Haidl, and Sergei Gorlatch. Parallelizing Branch-and-Bound on GPUs for Optimization of Multiproduct Batch Plants. In Victor Malyshkin, editor, Parallel Computing Technologies (PaCT 2015), volume 9251 of Lecture Notes in Computer Science, pages 324–337, Cham, 2015. Springer. doi: 10.1007/978-3-319-21909-7_33. [95] Andrey Borisenko, Michael Haidl, and Sergei Gorlatch. A GPU parallelization of branchand-bound for multiproduct batch plants optimization. The Journal of Supercomputing, 73 (2):639–651, 2017. doi: 10.1007/s11227-016-1784-x. [96] Jan Gmys, Mohand Mezmaz, Nouredine Melab, and Daniel Tuyttens. A GPU-based Branchand-Bound Algorithm using Integer-Vector-Matrix Data Structure. Parallel Computing, 59:119–139, 2016. doi: 10.1016/j.parco.2016.01.008. 92

[97] Nouredine Melab, Jan Gmys, Mohand Mezmaz, and Daniel Tuyttens. Multi-core versus Many-core Computing for Many-task Branch-and-Bound Applied to Big Optimization Problems. Future Generation Computer Systems, 82:472–481, 2018. [98] Ashkan Tousimojarad and Wim Vanderbauwhede. Number of Tasks, not Threads, is Key. In 2015 23rd Euromicro International Conference on Parallel, Distributed, and Network-Based Processing, pages 128–136. IEEE, 2015. doi: 10.1109/PDP.2015.81. [99] A. B. Borisenko and S. V. Karpushkin. Using OpenMP for the Optimal Equipment Selection of Multiproduct Plants. Vychislitel’nye Tekhnologii [Computational Technologies], 20(4):17–28, 2015. In Russian. [100] OpenMP Architecture Review Board. The OpenMP API Specification for Parallel Programming, 2022. URL https://www.openmp.org/. Accessed 2026-10-18. [101] K. Aida, W. Natsume, and Y. Futakata. Distributed computing with hierarchical masterworker paradigm for parallel branch and bound algorithm. In CCGrid 2003. 3rd IEEE/ACM International Symposium on Cluster Computing and the Grid, 2003. Proceedings., pages 156–163. IEEE, 2003. ISBN 0-7695-1919-9. doi: 10.1109/CCGRID.2003.1199364. [102] Hinde Lilia Bouziane, Christian Pérez, and Thierry Priol. Extending Software Component Models with the Master-Worker Paradigm. Parallel Computing, 36(2-3):86–103, 2010. doi: 10.1016/j.parco.2009.12.012. [103] Message Passing Interface Forum. Official MPI (Message Passing Interface) Standards Documents, Errata, 2022. URL https://www.mpi-forum.org/. Accessed 2026-10-18. [104] Kento Aida, Yoshiaki Futakata, and Tomotaka Osumi. Parallel Branch and Bound Algorithm with the Hierarchical Master-Worker Paradigm on the Grid. IPSJ Digital Courier, 2:584–597, 2006. doi: 10.2197/ipsjdc.2.584. [105] Abel Castellanos, Andreu Moreno, Joan Sorribes, and Tomas Margalef. Performance Model for Master/Worker Hybrid Applications on Multicore Clusters. In 2013 IEEE 10th International Conference on High Performance Computing and Communications & 2013 IEEE International Conference on Embedded and Ubiquitous Computing, pages 210–217. IEEE, 2013. ISBN 978-0-7695-5088-6. doi: 10.1109/HPCC.and.EUC.2013.39. [106] Wen-mei W. Hwu. Heterogeneous System Architecture: A new compute platform infrastructure. Morgan Kaufmann, 2015. ISBN 978-0-12-800386-2. [107] A. B. Borisenko. Development and Study of Parallel Metaheuristic Algorithms on Heterogeneous Computing Systems. In Proceedings of the Scientific Seminar of the “Mikhail Lomonosov” and “Immanuel Kant” Programme Fellows 2016/2017, pages 57–62, Moscow, 2017. Flinta. In Russian. [108] John D. Owens, David Luebke, Naga Govindaraju, Mark Harris, Jens Krüger, Aaron E. Lefohn, and Timothy J. Purcell. A Survey of General-Purpose Computation on Graphics Hardware. Computer Graphics Forum, 26(1):80–113, 2007. doi: 10.1111/j.1467-8659.2007 .01012.x. [109] NVIDIA Corporation. CUDA Programming Guide, 2026. URL https://docs.nvidi a.com/cuda/cuda-programming-guide/pdf/cuda-programming-guide.pdf. Accessed 2026-10-18. [110] A. V. Boreskov and A. A. Kharlamov. Fundamentals of CUDA Technology. DMK Press, Moscow, 2010. ISBN 978-5-94074-578-5. In Russian. 93

[111] A. M. Kazennov. Basic Concepts of CUDA Technology. Computer Research and Modeling, 2(3):295–308, 2010. doi: 10.20537/2076-7633-2010-2-3-295-308. In Russian. [112] J. Sanders and E. Kandrot. CUDA by Example: An Introduction to General-Purpose GPU Programming. Addison-Wesley, Upper Saddle River, NJ, 2010. ISBN 978-0-13-138768-3. [113] Javier Diaz, Camelia Munoz-Caro, and Alfonso Nino. A Survey of Parallel Programming Models and Tools in the Multi and Many-Core Era. IEEE Transactions on Parallel and Distributed Systems, 23(8):1369–1386, 2012. doi: 10.1109/TPDS.2011.308. [114] David B. Kirk and Wen-mei W. Hwu. Programming massively parallel processors: a hands-on approach. Morgan Kaufmann, 2016. ISBN 978-0-12-3914187. [115] Inc The Khronos Group. The Open Standard for Parallel Programming of Heterogeneous Systems, 2026. URL https://www.khronos.org/opencl/. Accessed 2026-10-18. [116] Aaftab Munshi, Benedict Gaster, Timothy G. Mattson, James Fung, and Dan Ginsburg. OpenCL Programming Guide. Addison-Wesley Professional, 2011. [117] Jason Zink, Matt Pettineo, and Jack Hoxley. Practical rendering and computation with Direct3D 11. AK Peters/CRC Press, 2016. ISBN 978-1568817200. [118] Ronald L. Graham, Donald E. Knuth, and Oren Patashnik. Concrete Mathematics: A Foundation for Computer Science. Addison-Wesley, Reading, Massachusetts, second edition, 1994. ISBN 0-201-55802-5. [119] A. B. Borisenko and S. V. Karpushkin. A Parallel Algorithm for the Optimal Equipment Selection of Multiproduct Plants. Prikladnaya Informatika [Journal of Applied Informatics], (3(45)):56–68, 2013. In Russian. [120] NVIDIA Corporation. NVIDIA’s Next Generation CUDA Compute Architecture: Kepler GK110/210, 2014. URL https://www.nvidia.com/content/dam/en-zz/Solutions/Da ta-Center/tesla-product-literature/NVIDIA-Kepler-GK110-GK210-Architectur e-Whitepaper.pdf. Accessed 2026-10-18. [121] Pedro Bruel, Marcos Amarís, and Alfredo Goldman. Autotuning CUDA compiler parameters for heterogeneous applications using the OpenTuner framework. Concurrency and Computation: Practice and Experience, 29(22):e3973, 2017. doi: 10.1002/cpe.3973. [122] Yingchao Huang and Dong Li. Performance Modeling for Optimal Data Placement on GPU with Heterogeneous Memory Systems. In 2017 IEEE International Conference on Cluster Computing (CLUSTER), pages 166–177. IEEE, 2017. ISBN 978-1-5386-2326-8. doi: 10.1109/CLUSTER.2017.42. [123] Tianyi David Han and Tarek S. Abdelrahman. Reducing branch divergence in GPU programs. In Proceedings of the Fourth Workshop on General Purpose Processing on Graphics Processing Units, pages 1–8. ACM, 2011. doi: 10.1145/1964179.1964184. [124] Imen Chakroun, Mohand Mezmaz, Nouredine Melab, and Ahcène Bendjoudi. Reducing thread divergence in a GPU-accelerated branch-and-bound algorithm. Concurrency and Computation: Practice and Experience, 25(8):1121–1136, 2013. doi: 10.1002/cpe.2931. [125] Huanxin Lin, Cho-Li Wang, and Hongyuan Liu. On-GPU Thread-Data Remapping for Branch Divergence Reduction. ACM Transactions on Architecture and Code Optimization, 15(3):1–24, 2018. doi: 10.1145/3242089.

94

[126] Junji Fukuhara and Munehiro Takimoto. Branch Divergence Reduction Based on Code Motion. Journal of Information Processing, 28:302–309, 2020. doi: 10.2197/ipsjjip.28.302. [127] Tao Zhang, Wei Shu, and Min-You Wu. Optimization of N-Queens Solvers on Graphics Processors. In Advanced Parallel Processing Technologies (APPT 2011), volume 6965 of Lecture Notes in Computer Science, pages 142–156, Berlin, Heidelberg, 2011. SpringerVerlag. ISBN 9783642241505. doi: 10.1007/978-3-642-24151-2_11. [128] Kamil Rocki and Reiji Suda. Parallel Minimax Tree Searching on GPU. In Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), volume 6067, pages 449–456. Springer, 2010. ISBN 364214389X. doi: 10.1007/978-3-642-14390-8_47. [129] T. S. Lee, Ghorbanali Moslemipour, T. O. Ting, and Dirk Rilling. A Novel Hybrid ACO/SA Approach to Solve Stochastic Dynamic Facility Layout Problem (SDFLP). In Emerging Intelligent Computing Technology and Applications, pages 100–108. Springer, 2012. doi: 10.1007/978-3-642-31837-5_15. [130] Christine Solnon. Ant Colony Optimization and Constraint Programming. John Wiley & Sons, Inc., Hoboken, NJ USA, 2013. ISBN 9781118557563. doi: 10.1002/9781118557563. [131] Francesca Rossi, Peter Van Beek, and Toby Walsh. Handbook of constraint programming. Elsevier, 2006. ISBN 9780444527264. [132] O. A. Shcherbina. Metaheuristic Algorithms for Combinatorial Optimization Problems (a Review). Tavricheskii Vestnik Informatiki i Matematiki [Taurida Journal of Computer Science Theory and Mathematics], (1(24)):56–72, 2014. In Russian. [133] A. N. Shabelnikov. Probabilistic Algorithms for Solving Dispatching Problems at Marshalling Yards and Humps. Izvestiya SFedU. Engineering Sciences, (7(168)):250–262, 2015. In Russian. [134] David Hibler. The GA-ACO Method Applied to Engineering Design. Procedia Computer Science, 12:152–157, 2012. doi: 10.1016/j.procs.2012.09.046. [135] A. V. Panteleev, D. V. Metlitskaya, and E. A. Aleshina. Global Optimization Methods. Metaheuristic Strategies and Algorithms. Vuzovskaya Kniga, Moscow, 2013. ISBN 978-59502-0743-3. In Russian. [136] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by Simulated Annealing. Science, 220(4598):671–680, 1983. doi: 10.1126/science.220.4598.671. [137] A. S. Tikhomirov. On Fast Variants of the Annealing Algorithm. Stokhasticheskaya Optimizatsiya v Informatike [Stochastic Optimization in Informatics], 5:65–90, 2009. In Russian. [138] Jayaraman Valadi and Patrick Siarry, editors. Applications of Metaheuristics in Process Engineering. Springer International Publishing, Cham, 2014. ISBN 978-3-319-06507-6. doi: 10.1007/978-3-319-06508-3. [139] Alexander G. Nikolaev and Sheldon H. Jacobson. Simulated Annealing. In Handbook of Metaheuristics, pages 1–39. Springer, 2010. doi: 10.1007/978-1-4419-1665-5_1. [140] S. S. Skiena. The Algorithm Design Manual. Springer, London, 2nd edition, 2008. doi: 10.1007/978-1-84800-070-4.

95

[141] O. A. Shcherbina. Constraint Satisfaction and Constraint Programming. Intellektual’nye Sistemy [Intelligent Systems], 15(1–4):53–170, 2011. In Russian. [142] Marco Dorigo and Thomas Stützle. Ant Colony Optimization: Overview and Recent Advances. In Handbook of metaheuristics, pages 311–351. Springer, 2019. doi: 10.1007/97 8-3-319-91086-4_10. [143] S. D. Shtovba. Ant Algorithms. Exponenta Pro. Matematika v Prilozheniyakh [Exponenta Pro. Mathematics in Applications], (4):70–75, 2003. In Russian. [144] Yu. Zaychenko and N. Murga. A Study of Ant Colony Optimization Algorithms for the Travelling Salesman Problem. Information Models and Analyses, 2(4):370–384, 2013. In Russian. [145] Kavishi Agarwal, Akshita Sinha, and M. Hima Bindu. A Novel Hybrid Approach to N-Queen Problem. In Advances in Computer Science, Engineering & Applications, pages 519–527. Springer, 2012. doi: 10.1007/978-3-642-30157-5_52. [146] Salabat Khan, Mohsin Bilal, M. Sharif, Malik Sajid, and Rauf Baig. Solution of n-Queen problem using ACO. In 2009 IEEE 13th International Multitopic Conference, pages 1–5. IEEE, 2009. ISBN 978-1-4244-4872-2. doi: 10.1109/INMIC.2009.5383157. [147] Tony White, Simon Kaegi, and Terri Oda. Revisiting Elitism in Ant Colony Optimization. In Genetic and Evolutionary Computation – GECCO 2003, pages 122–133. Springer, 2003. doi: 10.1007/3-540-45105-6_11. [148] Thomas Stützle, Manuel López-Ibáñez, Paola Pellegrini, Michael Maur, Marco Montes de Oca, Mauro Birattari, and Marco Dorigo. Parameter Adaptation in Ant Colony Optimization. In Youssef Hamadi, Eric Monfroy, and Frédéric Saubion, editors, Autonomous Search, pages 191–215. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011. ISBN 978-3642-21434-9. doi: 10.1007/978-3-642-21434-9_8. [149] Audrey Delévacq, Pierre Delisle, Marc Gravel, and Michaël Krajecki. Parallel Ant Colony Optimization on Graphics Processing Units. Journal of Parallel and Distributed Computing, 73(1):52–61, 2013. doi: 10.1016/j.jpdc.2012.01.003. [150] Anna Syberfeldt and Tom Ekblom. A Comparative Evaluation of the GPU vs The CPU for Parallelization of Evolutionary Algorithms Through Multiple Independent Runs. International Journal of Computer Science and Information Technology, 9(3):01–14, 2017. doi: 10.5121/ijcsit.2017.9301. [151] Fariha Nosheen, Sadia Bibi, and Salabat Khan. Ant Colony Optimization based scheduling algorithm. In 2013 International Conference on Open Source Systems and Technologies, pages 18–22. IEEE, 2013. ISBN 978-1-4799-2046-4. doi: 10.1109/ICOSST.2013.6720599. [152] Arushi Gupta and Smriti Srivastava. Comparative Analysis of Ant Colony and Particle Swarm Optimization Algorithms for Distance Optimization. Procedia Computer Science, 173:245–253, 2020. doi: 10.1016/j.procs.2020.06.029. [153] NVIDIA Corporation. Random Number Generation on NVIDIA GPUs (cuRAND), 2026. URL https://developer.nvidia.com/curand. Accessed 2026-10-18. [154] Vasiliy N. Leonenko, Nikolai V. Pertsev, and Marc Artzrouni. Using High Performance Algorithms for the Hybrid Simulation of Disease Dynamics on CPU and GPU. Procedia Computer Science, 51:150–159, 2015. doi: 10.1016/j.procs.2015.05.214.

96

[155] Tair Askar, Bekdaulet Shukirgaliyev, Martin Lukac, and Ernazar Abdikamalov. Evaluation of Pseudo-Random Number Generation on GPU Cards. Computation, 9(12):142, 2021. doi: 10.3390/computation9120142. [156] Md. Maruf Hussain, Hiroshi Hattori, and Noriyuki Fujimoto. A CUDA Implementation of the Standard Particle Swarm Optimization. In 2016 18th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC), pages 219–226. IEEE, 2016. ISBN 978-1-5090-5707-8. doi: 10.1109/SYNASC.2016.043. [157] Andrey Borisenko and Sergei Gorlatch. Parallelizing Metaheuristics for Optimal Design of Multiproduct Batch Plants on GPU. In Victor Malyshkin, editor, Parallel Computing Technologies (PaCT 2017), volume 10421 of Lecture Notes in Computer Science, pages 405–417, Cham, 2017. Springer. doi: 10.1007/978-3-319-62932-2_39. [158] Andrey Borisenko and Sergei Gorlatch. Comparing GPU-parallelized metaheuristics to branch-and-bound for batch plants optimization. The Journal of Supercomputing, 75(12): 7921–7933, 2019. doi: 10.1007/s11227-018-2472-9. [159] Paola Pellegrini, Thomas Stützle, and Mauro Birattari. A Critical Analysis of Parameter Adaptation in Ant Colony Optimization. Swarm Intelligence, 6(1):23–48, 2012. doi: 10.1007/s11721-011-0061-0. [160] Yasemin Eryoldaş and Alptekin Durmuşoglu. A Literature Survey on Offline Automatic Algorithm Configuration. Applied Sciences, 12(13):6316, 2022. doi: 10.3390/app12136316. [161] Mehdi Fallahi, Somayeh Amiri, and Masoud Yaghini. A parameter tuning methodology for metaheuristics based on design of experiments. International Journal of Engineering and Technology Sciences, 2(6):497–521, 2014. [162] Jianyong Sun, Qingfu Zhang, and Xin Yao. Meta-Heuristic Combining Prior Online and Offline Information for the Quadratic Assignment Problem. IEEE Transactions on Cybernetics, 44(3):429–444, 2014. doi: 10.1109/TCYB.2013.2256892. [163] Ren Diao and Qiang Shen. Deterministic parameter control in Harmony Search. In 2010 UK Workshop on Computational Intelligence (UKCI), pages 1–7. IEEE, 2010. ISBN 978-1-4244-8774-5. doi: 10.1109/UKCI.2010.5625576. [164] Michael Maur, Manuel Lopez-Ibanez, and Thomas Stutzle. Pre-scheduled and adaptive parameter variation in MAX-MIN Ant System. In IEEE Congress on Evolutionary Computation, pages 1–8. IEEE, 2010. ISBN 978-1-4244-6909-3. doi: 10.1109/CEC.2010.5 586332. [165] Efrén Mezura-Montes and Ana Gabriela Palomeque-Ortiz. Self-adaptive and Deterministic Parameter Control in Differential Evolution for Constrained Optimization. In ConstraintHandling in Evolutionary Optimization, pages 95–120. Springer, 2009. doi: 10.1007/978-3 -642-00619-7_5. [166] Eduardo De Mendonca Mesquita, Renato Coral Sampaio, Helon Vicente Hultmann Ayala, and Carlos Humberto Llanos. Recent Meta-Heuristics Improved by Self-Adaptation Applied to Nonlinear Model-Based Predictive Control. IEEE Access, 8:118841–118852, 2020. doi: 10.1109/ACCESS.2020.3005318. [167] Vincent A. Cicirello. Self-Tuning Lam Annealing: Learning Hyperparameters While Problem Solving. Applied Sciences, 11(21):9828, 2021. doi: 10.3390/app11219828. [168] André S. Santos, Ana M. Madureira, and Leonilde R. Varela. A Self-Parametrization Framework for Meta-Heuristics. Mathematics, 10(3):475, 2022. doi: 10.3390/math10030475. 97

[169] Davide Anghinolfi, Antonio Boccalatte, Massimo Paolucci, and Christian Vecchiola. Performance Evaluation of an Adaptive Ant Colony Optimization Applied to Single Machine Scheduling. In Asia-Pacific Conference on Simulated Evolution and Learning, volume 5361 of Lecture Notes in Computer Science, pages 411–420. Springer, 2008. doi: 10.1007/978-3-540-89694-4_42. [170] Felix Dobslaw. Recent development in automatic parameter tuning for metaheuristics. In Proceedings of the 19th Annual Conference of Doctoral Students-WDS 2010, pages 54–63, 2010. [171] Jerry Swan, Steven Adriaensen, Alexander E. I. Brownlee, Kevin Hammond, Colin G. Johnson, Ahmed Kheiri, Faustyna Krawiec, J. J. Merelo, Leandro L. Minku, Ender Özcan, Gisele L. Pappa, Pablo García-Sánchez, Kenneth Sörensen, Stefan Voß, Markus Wagner, and David R. White. Metaheuristics ”In the Large”. European Journal of Operational Research, 297(2):393–406, 2022. doi: 10.1016/j.ejor.2021.05.042. [172] Felix Dobslaw. A parameter tuning framework for metaheuristics based on design of experiments and artificial neural networks. In International Conference on Computer Mathematics and Natural Computing, volume 64, pages 213–216. WASET, 2010. [173] Arif Arin, Ghaith Rabadi, and Resit Unal. Comparative studies on design of experiments for tuning parameters in a genetic algorithm for a scheduling problem. International Journal of Experimental Design and Process Optimisation, 2(2):102–124, 2011. doi: 10.150 4/IJEDPO.2011.040262. [174] Mauro Birattari, Thomas Stützle, Luis Paquete, and Klaus Varrentrapp. A Racing Algorithm for Configuring Metaheuristics. In GECCO 2002: Proceedings of the Genetic and Evolutionary Computation Conference, pages 11–18, San Francisco, 2002. Morgan Kaufmann. [175] Manuel López-Ibáñez, Jérémie Dubois-Lacoste, Leslie Pérez Cáceres, Mauro Birattari, and Thomas Stützle. The irace package: Iterated racing for automatic algorithm configuration. Operations Research Perspectives, 3:43–58, 2016. doi: 10.1016/j.orp.2016.09.002. [176] Frank Hutter, Holger H. Hoos, Kevin Leyton-Brown, and Thomas Stützle. ParamILS: An Automatic Algorithm Configuration Framework. Journal of Artificial Intelligence Research, 36:267–306, 2009. [177] Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown. Sequential Model-Based Optimization for General Algorithm Configuration. In Learning and Intelligent Optimization (LION 5), volume 6683 of Lecture Notes in Computer Science, pages 507–523. Springer, 2011. doi: 10.1007/978-3-642-25566-3_40. [178] Áthila R. Trindade and Felipe Campelo. Tuning metaheuristics by sequential optimisation of regression models. Applied Soft Computing, 85:105829, 2019. doi: 10.1016/j.asoc.2019. 105829. [179] Peng Li and Hua Zhu. Parameter Selection for Ant Colony Algorithm Based on Bacterial Foraging Algorithm. Mathematical Problems in Engineering, 2016:1–12, 2016. doi: 10.115 5/2016/6469721. [180] Vladimir N. Temlyakov and Pavel Zheltov. On performance of greedy algorithms. Journal of Approximation Theory, 163(9):1134–1145, 2011. doi: 10.1016/j.jat.2011.03.009. [181] Urfat Nuriyev, Onur Ugurlu, and Fidan Nuriyeva. Self-Organizing Iterative Algorithm for Travelling Salesman Problem. IFAC-PapersOnLine, 51(30):268–270, 2018. doi: 10.1016/j. ifacol.2018.11.299. 98

[182] Gözde Kizilateş and Fidan Nuriyeva. On the Nearest Neighbor Algorithms for the Traveling Salesman Problem. In Advances in Computational Science, Engineering and Information Technology, pages 111–118. Springer, 2013. doi: 10.1007/978-3-319-00951-3_11. [183] Frumen Olivas, Fevrier Valdez, and Oscar Castillo. Dynamic parameter adaptation in Ant Colony Optimization using a fuzzy system for TSP problems. In IFSA-EUSFLAT, pages 765–770, 2015. doi: 10.2991/ifsa-eusflat-15.2015.108. [184] Oscar Castillo, Héctor Neyoy, José Soria, Patricia Melin, and Fevrier Valdez. A new approach for dynamic fuzzy logic parameter tuning in Ant Colony Optimization and its application in fuzzy control of a mobile robot. Applied Soft Computing, 28:150–159, 2015. doi: 10.1016/j.asoc.2014.12.002. [185] Eduardo Batista de Moraes Barbosa and Edson Luiz França Senne. Improving the FineTuning of Metaheuristics: An Approach Combining Design of Experiments and Racing Algorithms. Journal of Optimization, 2017:1–7, 2017. doi: 10.1155/2017/8042436. [186] E. S. Skakov and V. N. Malysh. Parameter meta-optimization of metaheuristics of solving specific NP-hard facility location problem. Journal of Physics: Conference Series, 973: 012063, 2018. doi: 10.1088/1742-6596/973/1/012063. [187] Marius Lindauer, Holger H. Hoos, Kevin Leyton-Brown, and Torsten Schaub. Automatic construction of parallel portfolios via algorithm configuration. Artificial Intelligence, 244: 272–290, 2017. [188] Andrey Borisenko and Sergei Gorlatch. Optimizing a GPU-Parallelized Ant Colony Metaheuristic by Parameter Tuning. In Victor Malyshkin, editor, Parallel Computing Technologies (PaCT 2019), volume 11657 of Lecture Notes in Computer Science, pages 151–165, Cham, 2019. Springer. doi: 10.1007/978-3-030-25636-4_12. [189] Andrey Borisenko and Sergei Gorlatch. Efficient GPU-parallelization of batch plants design using metaheuristics with parameter tuning. Journal of Parallel and Distributed Computing, 154:74–81, 2021. doi: 10.1016/j.jpdc.2021.03.012. [190] Carla P. Gomes, Bart Selman, Nuno Crato, and Henry Kautz. Heavy-Tailed Phenomena in Satisfiability and Constraint Satisfaction Problems. Journal of Automated Reasoning, 24(1–2):67–100, 2000. [191] Enrique Alba. Parallel Evolutionary Algorithms Can Achieve Super-Linear Performance. Information Processing Letters, 82(1):7–13, 2002.

99

Related documents

• Multi-fidelity batch Bayesian optimization for bioprocess development across scales#996871
• OverlapAdvisor: Modeling and Adapting Computation-Communication Overlap for Distributed MoE Serving#1636358
• OverlapAdvisor: Modeling and Adapting Computation-Communication Overlap for Distributed MoE Serving#1605203
• A Guide to Bayesian Optimization in Bioprocess Engineering#1278781
• Multi-objective production planning of new and remanufactured products in hybrid production system#643342