A Practical Methodology for Reproducible Experimentation: An Application to the Double-Row Facility Layout Problem

There exists a growing concern in the research community about the lack of methodological frameworks to compare stochastic algorithms. The list of problems ranges from the election of instances (Bartz-Beielstein et al., 2020) to ensuring the experimental reproducibility (Kendall et al., 2016). Evolutionary computation (Bäck et al., 1997) in particular, and metaheuristics (Glover and Kochenberger, 2006) in general, are among the most renowned stochastic algorithms. In this kind of procedure, results are supported by empirical studies. Therefore, the ability to reach similar results (both in quality and computing time) by repeating an experiment performed by other researchers is the only way in which the research community can reach a consensus about an empirical finding. Although a set of guidelines for optimization studies in software engineering exists as part of the ACM SGSOFT empirical standards effort (Ralph et al., 2020), many scientific works are not reproducible at all (Baker, 2016). Moreover, reproducibility issues do not only affect researchers when reproducing experiments conducted by other researchers, but also when reproducing their own previous experiments.

In this article, we mainly focus on reproducibility in the context of stochastic procedures. It is worth mentioning that there does not exist a well established terminology in this matter. Indeed, terms such as reproducibility and replicability are used almost synonymously (López-Ibáñez et al., 2021). We refer the reader to Plesser (2018) for further details about the evolution of the terminology in this area.

The first studies of reproducibility can be traced back to the 1990s (Moscato and Norman, 1992). In the same line, Gent et al. (1997) claimed that the empirical study of algorithms is a relatively “immature” field, proposing (colloquially) some “thing to do” and “thing to not do” in order to ensure the experimental reproducibility. Later, Johnson (2002) provided a useful guide to computer scientists about how experimental analysis of algorithms can be performed and, then, described. A similar analysis was conducted in Eiben and Jelasity (2002) but focused on evolutionary algorithms. The Association of Computing Machinery (ACM, 2021) has recently specified different degrees of reproducibility by considering the concepts of artifact and measurement. Specifically, an artifact is defined as a digital object that was either created by researchers to be used as part of the study or generated by the experiment itself. Examples of artifacts in the context of stochastic procedures would be source/executable code of an algorithm or data/code required to fully specify benchmark instances. A measurement is introduced as an analogy to physical experiments. In the context of stochastic procedures, it represents the raw data (objective function values, runtimes, etc.) that results from an experiment. Notice that depending on the level of abstraction being studied, sometimes reporting aggregated results is more adequate; for example, summary statistics such as means and standard errors, computational effort evaluated as cycles, function evaluations, or iteration counters.

Based on these definitions, ACM introduced three new concepts: repeatability, reproducibility, and replicability (ACM, 2021). However, more recently, a redefinition of this classification was proposed in López-Ibáñez et al. (2021), also adding the concept of generalizability. This new classification is based on the definition of fixed and random factors, which are determined by the researcher and given by elements like random seeds, respectively. The Turing Way project (Arnold et al., 2019) introduced a complementary approach. Specifically, it considers reproducibility (same analysis performed on the same dataset consistently produces the same answer), replicability (same analysis performed on different datasets produces qualitatively similar answers), robustness (different analysis applied to the same dataset produces a qualitatively similar answer to the same question), and generalizability (the combination of replicability and robustness). This classification has been mainly adopted in the machine learning community (Pineau et al., 2021). An alternative approach is introduced in Stodden et al. (2014), which relates to the availability of the code, data, and all the details of the implementation and experimental setup that allow obtaining the published result. It also considers statistical reproducibility, which is concerned with validating the results of repeated experiments by means of statistical assessments.

There exists a vast literature (see López-Ibáñez et al., 2021, for a recent survey) asking for general guidelines that aim to ensure, assess, and encourage the reproducibility of experiments; for instance, by publishing permanently accessible, complete, and useful artifacts. Unfortunately, this is particularly difficult in stochastic procedures, since one must deal not only with differences on hardware and/or software for reproducing experiments, but also with algorithm randomness. Replicating the exact conditions of previously reported experiments might be impossible, even for those researchers who originally proposed the corresponding experiment. This situation arises when, for example, the original hardware is not available anymore, the version of some software libraries is not specified, or some sources of randomness are not repeatable. In that case, repeatability and reproducibility are unreachable. Replicability may still be achievable, but unfortunately, it is not an easy task since replicability requires high-level descriptions of artifacts with enough detail to enable their independent implementation and a careful selection of measurements with their stated precision and confidence levels. Only if those prerequisites are fulfilled, can other researchers unequivocally conclude if the replication of an experiment is confirmed.

In this article, we propose a methodological approach to increase the reproducibility of empirical results. It is based on an automated characterization and selection of test instances and an automated configuration of the proposed algorithm. In addition, source code, instances, and executable artifacts are publicly available in Zenodo. With this main contribution, we validate the proposed methodology with a well-known combinatorial optimization problem. Specifically, we propose a new metaheuristic method for the DRFLP, which is able to obtain better results than the state-of-the-art algorithms using shorter computation times. In particular, the proposal obtains all best-known values in a fraction of the time required by previous methods. The new algorithm is not an outcome of the methodology, but the reproducibility of the conducted experimentation is.

The rest of the article is organized as follows. The proposed methodology for reproducibility is described in Section 2. The definition of the target problem and the algorithmic proposals are detailed in Sections 3 and 4, respectively. The experimental results and the application of the proposed methodology to the target problem are explained in Section 5, and a thorough description of the provided artifacts is given in Section 6. Finally, conclusions and future work are drawn in Section 7.

The main goal of this article is to encourage reproducibility efforts for optimization problems being solved with metaheuristics. In order to do so, we propose an empirical methodology based on minimizing the number of decisions taken by researchers (i.e., automatizing the majority of them). It is worth mentioning that this methodology is suitable not only for any evolutionary algorithm, but also for metaheuristic procedures in general.

The proposed methodology is in line with the one described in López-Ibáñez et al. (2021). As principal outcomes, it allows researchers to conduct studies about repeatability, reproducibility, replicability, and generalizability. For the sake of completeness, we include the definition of these concepts, as well as the definition of random and fixed factors in the context of performance assessment of stochastic algorithms:

• Random factors correspond to those elements whose value may belong to a certain distribution, but only a random subset of values will be evaluated in the experiment. The typical random factors are the random seeds.

• Fixed factors, as in evolutionary computation, include target instances, parameter values, time limits, and so on. Converting fixed factors into random factors allows the generalization of conclusions from an experiment.

• Repeatability consists in repeating the proposed experiment, obtaining the very same results. This kind of study requires the availability of the same instances, artifacts, and run environment as the original study. As will be later described in more detail, the software framework used is able to generate a Docker container with the full run environment, including the same random factors, in order to guarantee that the results are identical after each execution.

• Reproducibility analyzes the influence of random factors over the results of an experiment. In the proposed methodology, an independent artifact is provided which automates the benchmark instance selection. In addition, an external tool is used to select the parameter values of the algorithm. Finally, the algorithm and source code are provided as artifacts. Therefore, a researcher may vary the random seeds in any of these elements, including the benchmark selection process.

• Replicability compares the results from a new implementation of the artifacts using new random factors. In this case, since our methodology thoroughly describes the benchmark instances selection process and the algorithm configuration, and the source code is also available, any replication is possible.

• Generalizability comes after replicability studies, also analyzing new fixed factors. Again, since all the steps of the methodology are detailed, a researcher may change the fixed factors and compare the obtained results. Notice that elements of the benchmark selection like the features of the instances are also fixed factors which can be changed.

Next, each one of the proposed steps of our methodology is thoroughly described, as well as the selected implementation for the target problem.

Computer scientists have mainly considered two different strategies to deal with instances, which are based on separating them into two different subsets. The first one divides the set of instances between “training set” (used to configure/tune the proposed algorithm) and “testing set” (used to compare the proposed algorithm with the state-of-the-art procedures), the intersection of these two sets being empty. This strategy is the one mainly adopted by the machine learning (ML) community. The second strategy also divides the set of instances into two sets. On the one hand, there is the preliminary instance set or benchmark instances devoted to tuning/configuring the algorithm (equivalent to the training set). However, the comparison of the proposed algorithm with the state-of-the-art procedures is conducted over the whole set of available instances. This strategy is the one followed by the heuristic optimization community. As aforementioned, we validate the proposed empirical methodology by proposing a metaheuristic method to solve a specific optimization problem. Therefore, we consider an approach commonly used in the heuristic optimization community.

As a departing hypothesis, we consider that the larger the studied set of instances, the more precise the expected outcomes of the algorithm. Therefore, we presume that researchers will provide as many instances as possible for the corresponding experiments. According to the number of instances, it is customary to choose a representative percentage of the total number of instances. The design of this benchmark instance set is critical, affecting directly to the performance of an algorithm. As stated in De Souza et al. (2021), a common mistake in the process of tuning an algorithm occurs when the instances selected in the preliminary set do not represent the whole set of available instances. Therefore, the configurations of the algorithm may be specialized on “known instances” (those included in the preliminary set) and perform poorly on “unknown instances” (those not included in the preliminary set).

The size of the benchmark instance set is usually selected arbitrarily, and it depends, among other things, on the total number of instances of the problem, the time needed to perform the preliminary experiments, or the instances characteristics. Therefore, in this research, we do not go deeper into the decision of this value and let the researcher decide how many instances have to be selected in the benchmark set. However, we do address the decision of defining the benchmark instance set, that is, the choice of the instances that constitute the benchmark set.

We propose a semiautomatic benchmark selection process based on the four aforementioned principles: repeatability, reproducibility, replicability, and generalizability. To this end, we try to minimize the number of decisions based on the researcher's experience, by favoring those based on rational criteria. In particular, we adapt to the instance selection context some techniques commonly related to the statistical ML field. The idea of applying ML techniques in the metaheuristic domain to design efficient, effective, and robust algorithms has become increasingly popular over the last decades (Birattari et al., 2006; Olvera-López et al., 2010; Talbi, 2021). Among ML techniques, we propose the use of principal component analysis (PCA) and k-means clustering algorithms to select the benchmark instances.

As seen in Figure 1, the benchmark selection process receives as input the whole set of available instances of the problem. Additionally, it is provided a user-defined configuration that includes the number of instances to be selected, and the parameters required for the initialization of the PCA and k-means algorithms (seed, confidence level, number of repetitions, etc.). Then, four steps are sequentially executed: structural features characterization, principal component analysis, instance clusters generation, and benchmark instance ranking. These steps are explained below.

Structural Features Characterization. In this first step, the instances are read and processed to obtain, for each one of them, a set of metrics and features. This step needs the intervention of the researcher because it requires identifying and describing the features that are relevant for a given problem. Then, those features are extracted by the code implementation of the researcher. In many cases, understanding the problem domain helps to detect which features can be useful.

Principal Component Analysis. The previous step generates an output matrix of metrics and features, that is, a multivariate data set. However, a high number of instance features hinders the analysis of the data for extracting conclusions. This issue specifically gets worse when features have different scales or measurement units. Occasionally, this is not an issue at all, depending on the target problem. Therefore, we propose two well-known techniques to deal with it: standardization of the data and PCA.

Data standardization is a common requirement for ML algorithms and consists of transforming the data by removing the mean value of each feature and then scaling it by dividing each value by its standard deviation. Therefore, the application of data standardization results in a new set of values where the mean is zero and the standard deviation is one.

Once the data have been standardized, we perform the PCA. This technique is generally used to reduce the redundant information of a given set of variables. This procedure generates a new set of features, denoted as the principal components set, which retains much of the information from the original set, although the number of variables included has been significantly reduced (Jolliffe, 2005).

To determine the number of principal components, we use the explained variance ratio, that is, the percentage of variance that is attributed to each one of the selected components. Generally, we would like to select a number of principal components that explain around 90% of the data (Jolliffe and Cadima, 2016). We have implemented this process to run automatically. However, the percentage of attributed variance has to be established by the researcher in the input configuration file of the procedure.

Instance Clusters Generation. The result of the second step is a reduced set of variables that represent the features of the instances. In this third step, we use the k-means algorithm to classify the instances in clusters, according to the principal components previously obtained. The k-means algorithm is one of the most widely used clustering methods. It aims to partition a set of observations into k clusters, where more similar elements are classified in the same cluster (MacQueen et al., 1967).

The number of clusters is usually specified by the researcher. However, a simple and popular solution to determine the optimal number of clusters is the elbow method. The elbow method consists of executing the k-means algorithm for a determined range of k values. Then, for each cluster obtained in each one of the k executions, the average of the distortion score (the sum of the squared distances from each point to its assigned cluster center) is computed. Finally, the optimal value of k is determined visually as the point where the slope of the curve has a significant change (like an elbow) when plotting the obtained scores. In order to automate the process of determining the number of clusters and to avoid possible mistakes by the researcher, we have implemented the method proposed by Satopaa et al. (2011), which detects the optimal number of clusters and, thus, automating this process. Briefly, the authors proposed a formal definition for identifying a knee in discrete data sets based on the mathematical definition of curvature for continuous functions. We omit the details for the sake of brevity.

As a result of this third step, we obtain a classification of the instances into an appropriate set of clusters based on the initial structural features.

Benchmark Instances Ranking. The fourth and final step of the proposed methodology aims to determine the instances that will become part of the benchmark set. For this purpose, the k clusters are sorted in descending order according to their size, that is, the number of instances. Then, they are sequentially traversed, selecting from each of them the instance that minimizes the distance to its centroid. This automatic process is repeated until the desired number of instances of the benchmark set is reached. It is worth mentioning that the selected instances are removed within their cluster, so no instances can be repeated in the benchmark set, and if a cluster is empty, it is skipped and the next non-empty cluster is visited.

Hence, the selected benchmark instances represent the diversity of the complete set of instances according to their structural features. Moreover, the selection is not made by the researcher, but given by these automatic procedures, supported by the described methods.

The proposed methodology in this section has been implemented with Jupyter Notebooks using Python. Section 5.1 presents the application of the proposed methodology in the context of the DRFLP. All artifacts and source code are publicly available, as will be detailed in Section 6.

Generally, metaheuristic approaches require tuning several parameters. These parameters make the proposed algorithm robust, but also flexible, and have a great influence on its effectiveness and efficiency. The problem of finding the optimal parameters for an algorithm could be considered itself as a hard continuous optimization problem, where the objective is to optimize the performance of a given optimization algorithm over a set of preliminary instances.

Researchers have generally classified the tuning strategies in two types: offline parameter tuning and online parameter tuning. The first strategy consists of fixing the algorithm parameters before execution, while in the second strategy, the parameters are dynamically adapted during execution. This work is mainly focused on the first strategy, the offline parameter tuning.

Customarily, parameter tuning is done manually, and it includes a difficult and time-consuming test approach guided by the experience of the researcher. Particularly, manual tuning involves a preliminary exploration to find a suitable range of parameters. Then, tests are executed to evaluate only a single parameter at a time, while the others parameters are set to a rational value. The best value of the parameter tested is set for the next experiment. These systematic tests are performed until all algorithm parameters have been tuned (Crainic et al., 1993; Gendreau, 2003).

Although manual testing could be effective for adjusting a few parameters, it has many drawbacks when the algorithm to be optimized is complex. For instance, it has difficulties in adjusting algorithms in which there are cross effects between parameters, the design alternatives that are explored are limited by the number of experiments, the tuning process is irreproducible, or some configurations are discarded due to bad performance in a particular experiment, among others (Stützle and López-Ibáñez, 2019).

To mitigate the disadvantages of manual tuning, researchers have proposed advanced automatic tuning procedures that do not require their intervention. These procedures, usually denoted as automatic algorithm configuration (AAC) procedures, follow a common general scheme: first, they receive as input the parameters to configure of the algorithm, as well as their domain, range, and restrictions; then, they generate a list of candidate configurations of the algorithm which are evaluated over a benchmark set of instances; finally, the best configurations are returned. Examples of these AAC procedures are ParamILS (Hutter et al., 2009), sequential model-based algorithm configuration (SMAC) (Hutter et al., 2011), and irace (López-Ibáñez et al., 2016).

In this work, we propose the use of irace for automatically tuning the algorithm proposals. irace is an iterated racing method based on the Friedman-Race (F-Race) and the iterated F-Race (I/F-Race). Particularly, the irace program requires a user-defined input configuration that defines the behavior of the algorithm to be tuned. The evaluation of the configurations is done through a racing procedure where each configuration is evaluated for each instance of the benchmark set. The irace software has been implemented in R and is freely available at http://iridia.ulb.ac.be/irace/.

Section 5.2 describes the use of irace on the DRFLP, as well as the obtained results.

There exist several technical solutions that allow the portability of software to different environments such as virtual machines and containers like Docker, and platforms like Open Science Foundation or Code Ocean (Clyburne-Sherin et al., 2019). These technical solutions include everything needed to replicate an experiment (code, system tools, system libraries, and settings, and independence of the underlying operating system). Nevertheless, once this validation step is done, the preservation of code and data is more useful in the long term than the availability of a reproducible experimental environment. This is due to the rapid obsolesce of software/hardware. Even if the original study becomes nonreproducible due to the obsolescence of its original artifacts, studying their code and data could help future replication and generalization efforts.

According to the definitions from literature, an artifact is a digital object that was part of the research process, either created or being an external tool used in this process (ACM, 2021). In addition, instances, result files, etc. are also considered artifacts.

Therefore, taking into account the previously mentioned four different types of reproducible studies according to López-Ibáñez et al. (2021), our methodology proposes the generation of a number of mandatory artifacts, which are listed below:

• Problem instances and results files. These include not only the instance files to be read by the algorithm but also, if possible, explanatory text files regarding technical questions like, for instance, file format. In addition, plain text files (or spreadsheet files) are included with the obtained results.

• Source code of the proposed methods. These files include proposed software for instance processing (if needed), and proposed algorithms and methods to generate the final results with the metrics detailed in the work.

• Executable artifacts. These include binary or executable versions on a portable environment for all the elements described in the previous item.

• Analysis artifacts. These include generated diagrams, R scripts, or any tool used to generate them, with the raw results used to generate any related figure or summary table.

Giving access to these artifacts allows the development of any kind of reproducibility study. As will be described later, some of these artifacts were automatically generated by the framework we have developed.

While the main proposal of the work is to provide a practical methodology for reproducibility, we also describe how this methodology has been applied to a particular problem from the family of facility layout problems (FLP), where we have also improved the results from the state of the art.

The FLP family comprises a number of different optimization problems that are devoted to determine the optimal position of a given set of facilities in a particular layout. This family of problems has many interesting applications in flexible manufacturing systems (FMS), where facilities have to be positioned on a particular layout where, usually, an automated machine carries material among facilities (Satheesh Kumar et al., 2008). Hence, the objective of the optimization problem is to minimize a cost function that is related to the material handling cost between facilities.

Several layouts can be found in the related literature. For instance, the simplest layout, where all facilities are consecutively located, is known as the single-row facility layout problem (SRFLP) (Simmons, 1969; Rubio-Sánchez et al., 2016). Some other layouts are the multiple-row layout (Herrán et al., 2021), the circular layout (Hungerländer et al., 2020), and the T-row layout (Dahlbeck, 2021), among others.

In this article, we have selected one of the most common layouts in the literature, which is the use of two rows to arrange the facilities allowing free space separation between them, known as the DRFLP. Notice that the variant of this problem where no empty space is allowed between facilities is known as the space-free double-row facility layout problem (SF-DRFLP) or corridor allocation problem (CAP) (Amaral, 2012).

The DRFLP was firstly tackled from a heuristic point of view in Chung and Tanchoco (2010). In this work, five heuristics were compared with a branch-and-cut algorithm solving an MIP model. One of the heuristics was finally used to generate an initial solution for the exact algorithm, and its behavior was compared in a set of 16 instances from 6 to 36 facilities. The DRFLP was also studied in Amaral (2013), proposing a new mathematical model smaller than the one proposed by Chung and Tanchoco. The experimental results, however, do not compare the same instances as in the previous paper. Besides, some years later, the mathematical model was enhanced in Secchin and Amaral (2019) improving the results of Amaral (2013) in instances up to 16 facilities. In Chae and Regan (2020), the authors propose another mathematical model obtaining competitive results for instances up to 16 facilities. Following his work on this problem, Amaral recently published a work improving the mathematical model once again (Amaral, 2021) reporting better results than Chae and Regan (2020) and, finally, a different work proposing a two-phase algorithm that merges a heuristic method with a linear programming routine (Amaral, 2020). This is the paper with the highest number of instances and, in addition, with the largest sizes of the state of the art.

As seen, many works have proposed exact algorithms for this problem, and some of them have included heuristic methods. However, reproducibility has not been promoted in any of them. One of the reasons for this situation could be that there are no algorithmic parameters to tune due to the nature of the proposed methods. Besides, the small size of the instances studied in the state of the art prevents the design of an automated selection of instances. In addition, regarding the heuristic methods, there is no publicly available code or artifacts, and no effort has been made to group the studied instances in the most recent papers.

The DRFLP is an $NP$-hard optimization problem that consists in finding an optimal arrangement of rectangular facilities in two rows.

A solution to this optimization problem is to find an optimal mapping $π$ that assigns the set of facilities $F$ to the corresponding layout with two rows, together with the exact horizontal location of the center of each facility $i$ in the layout (abscissa $xi$⁠). Specifically, given a facility $i∈F$⁠, $πk(p)=i$ indicates that $i$ is located at position $p$ of the permutation corresponding to row $k$ in the layout $π$⁠. Hence, the complete layout is defined by two permutations of facilities, one for each row, $π={π1,π2}$⁠. It trivially holds that the number of facilities of a feasible solution is $|π|=|π1|+|π2|=n$⁠. For the example shown in Figure 2, $π1={F2,F3,F4,F8}$ and $π2={F7,F5,F6,F9,F1}$⁠. In addition to the permutation, the location of centers $xi$ must be provided since free space can be found in a solution.

As stated before, a solution to the DRFLP is defined by two permutations of facilities, one for each row, $π={π1,π2}$⁠, together with the exact horizontal location of the center of each facility $i$ in the layout (abscissa $xi$⁠). The location of the centers is required since an optimal solution of the DRFLP could have a free space between consecutive facilities located in the same row, and the permutation is not able to store this information without the centers.

If the length of each facility is an integer number greater than or equal to 1, it can be demonstrated that the allowed free spaces between consecutive facilities in the optimal solution (if they exist) is a multiple of 0.5 units.¹ Hence, a modification in the solution representation can be made in order to represent the complete layout (free spaces included) in the permutation.

In particular, we incorporate “dummy facilities” of length 0.5, whose weight between them and any other facility (including other dummy facilities) is equal to zero. Hence, adding a number of these “dummy facilities” to the original instance $I$⁠, makes it possible to find the optimal solution of an instance that could need free spaces between facilities by only dealing with permutations. This way, centers $xi$ are not considered as decision variables in this process.

Previous works on the DRFLP have designed heuristic methods which provided good-quality initial solutions for exact methods, as in Amaral (2020). On the contrary, our proposal is based on applying a metaheuristic method to tackle this problem. In particular, an iterated greedy algorithm is designed where the constructive and improvement phases are based on the greedy randomized adaptive search procedure (GRASP), followed by a destructive method. In addition, the algorithm is embedded in a multistart loop which is responsible for introducing spaces between facilities by means of “dummy facilities” (see Section 3.3).

Next, the algorithmic proposals are described, ending the section with the complete structure of the final algorithm.

We have selected the iterated greedy (IG) algorithm (Jacobs and Brusco, 1995) as the main component of our algorithmic proposal. It is a very fast and simple method based on consecutive partial destruction and reconstruction of a candidate solution. Due to its simple design and efficient performance, it has been applied to many different problems (Ruiz and Stützle, 2007; Lozano et al., 2011; Ruiz et al., 2019; Quintana et al., 2021).

Algorithm 1 shows the pseudocode of our IG proposal, which receives the following input parameters: $I$⁠, which is the target instance; $α1$ and $α2$⁠, which balance the random bias of the constructive methods; $imax$⁠, which corresponds to the number of iterations of the IG inner loop; and $β$⁠, which determines the damage ratio in the destructive method. As seen in step 1 of the pseudocode, an initial solution $π$ is generated from scratch. This initial solution is improved by a local search procedure and stored as the current best solution $π★$ in step 2. Then, the inner loop begins with the partial destruction of solution $π★$ in step 4. The resulting solution, $π'$⁠, is reconstructed in step 5 and then improved again by a local search in step 6, obtaining a new solution $π'''$⁠. Finally, steps 7 to 8 represent the acceptance criterion of our IG proposal, which is a simplified version of the one used in the original IG proposal from Ruiz and Stützle (2007). In particular, a reconstructed solution is accepted for the next iteration only if its cost function value is better than the current best solution $π★$⁠. Therefore, $π★$ is updated accordingly and returned by the algorithm in step 9 after the corresponding number of iterations $imax$⁠.

In order to either create an initial solution or reconstruct an incoming solution, we have designed two different approaches based on GRASP, which was proposed by Feo and Resende (1989) as a new method where greediness and randomness are balanced in the construction of a solution.

Algorithm 2 shows the pseudocode of our first constructive proposal, called C1. This method builds a solution by iteratively adding new facilities to the partial layout, which is initially empty. The method starts by adding a randomly selected facility $i$ at the beginning of the first row of an empty layout $π1(1)$ in steps 1 and 2. Then, after creating a candidate list $CL$ with the rest of facilities in step 3, the method enters a loop until all the facilities in $CL$ are allocated in the final layout $π$⁠, returned in step 13. To this purpose, an extended candidate list $ECL$ is created in step 5. Notice that this list contains all the combinations $(i,k,p)$ of facilities $i∈CL$ and the positions $1≤p≤|πk|+1$ in each row $1≤k≤2$ of the partial solution $π$ where these facilities can be inserted. Then, a greedy function $g$ is defined for the selection of a combination $(i,k,p)∈ECL$⁠. Specifically, $g$ evaluates the contribution to the objective function when a facility $i∈CL$ is inserted at position $p$ of row $k$⁠, resulting in a new partial solution $π'$⁠, $g(I,π,i,k,p)=F(I,π')-F(I,π)$⁠. Once the minimum and maximum values of this function are computed in steps 6 and 7, respectively, a threshold value $th$ is calculated in step 8 to generate the reduced candidate list $RCL$ from $ECL$ in 9. This list is populated with all the elements in $ECL$ whose greedy function $g$ is below the previously calculated threshold value $th$⁠, which depends on $α$⁠. Finally, the loop ends by randomly selecting one element $(i,k,p)$ from $RCL$ in step 10, and inserting the selected facility $i$ at the corresponding position $(k,p)$ of the layout $π$ in step 11. Notice that the insertion of a facility $i$ in position $(k,p)$ increases the center (abscissa) of all the facilities in positions from $p$ to $|πk|$ row $k$ in $li$ units. Once the facility $i$ is included in the partial layout $π$⁠, it is removed from $CL$ in step 12 and a new iteration stars computing the new $ECL$ at step 5.

As seen, this proposal follows the customary GRASP-constructive approach, which performs a greedy selection followed by a random choice of a candidate. An alternative to this approach is to exchange the random and greedy stages to perform random selections and, then, to select the best facility in terms of the greedy function. Therefore, we propose a new constructive procedure called $C1'$ where the main changes are in the while-loop. Specifically, steps 6 to 10 in Algorithm 2 are substituted by first taking at random a percentage $α$ of elements from $ECL$ (i.e., the number of selected elements is $α·|ECL|$⁠) and, finally, selecting the best one of these elements according to the corresponding greedy function. We have omitted the pseudocode description of this method for the sake of space, since the source code will be provided.

In addition to these GRASP-based constructive strategies, we propose a baseline method, called $C0$⁠, devoted to constructing a completely random solution. This method is equivalent to executing either $C1$ or $C2$ with $α=1$⁠. The goal of including this method is to evaluate whether the use of information about the problem (e.g., $C1$ and $C2$⁠) improves the outcomes.

As it is shown in Algorithm 1, our IG implementation includes a local search procedure to improve the solution generated in the construction phase. Specifically, we have used the best improvement strategy to explore the neighborhoods generated by insert moves. Moreover, we have also implemented the efficient computation of the cost associated to the neighbor solutions based on swap moves. See Herrán et al. (2021) for a detailed explanation of insert moves and the efficient exploration strategy.

In this work, we have implemented two different destructive methods depending on how the initial solution is constructed. The first method, called D1, randomly removes a number of facilities (controlled by the parameter $β$⁠) of the current solution $π$⁠.

Because we will provide the source code of our proposals, we have omitted the pseudocode of these methods for the sake of space.

In this article we propose a multistart iterated greedy algorithm (MS-IG) which we have modified in order to increasingly include “dummy facilities” into the corresponding instance $I$ while the current solution can be improved. Algorithm 4 shows the pseudocode of our final proposal. The algorithm begins by setting to 0 the number of “dummy facilities” of the instance to be solved (⁠$d$⁠) and the number of “dummy facilities” of the instance associated to the best current solution found (⁠$d★$⁠) in steps 1 and 2, respectively. It also calculates in step 3 the threshold value $th$ used for the termination criteria of the while loop. The threshold is defined as the maximum increment of fake facilities since last improvement, or, in other words, how many fake facilities can be added to a solution without improving before stopping. This value depends on the current instance, and it will be one of the parameters of the algorithm to tune. Next, the IG method described in Algorithm 1 is applied to solve the original instance $I$ (without including any “dummy facility”) in step 4. Then, the algorithm enters a while-loop, which applies $tmax$ times at each iteration the IG method, but using a new instance $I'$ containing additional “dummy facilities” (step 9). This new instance $I'$ is updated from $I$ in step 8 according to the new number of “dummy facilities” calculated in step 7. The strategy for increasing fake facilities can be easily customized. In this case, we propose two different implementations: linearly increment fake facilities after each iteration by a constant or exponentially increasing them using a function or sequence such as Fibonacci. Again, this will be another parameter to tune. After solving the new problem in step 9, steps 10 to 13 update the best current solution found $π$⁠, instance $I$ and number of “dummy facilities” $d★$ if some improvement on the objective function cost value is achieved in step 10. The loop in steps 5 to 13 iterates while the difference between the number of “dummy facilities” of the instance $I'$ solved at each iteration and the number of “dummy facilities” of the instance $I$ associated to the best current solution is below the threshold value $th$ calculated in step 3.

Once again, since the source code is available in Zenodo, we omit the pseudocode of GetThreshold, IncreaseDummyFacilities and IncludeDummyFacilities methods.

Once the proposed methodology is described, the target problem defined, and the algorithmic proposal detailed, this section shows the particular application of the methodology to the DRFLP. Similarly, we will also show the results obtained by our approach in relation to the state of the art on the target problem. Finally, since we propose a new set of larger instances, we have been forced to reimplement the algorithmic proposal from the state of the art because no code or binaries were provided by the authors. We then describe this replicability study as well as the obtained results.

All the algorithmic components have been implemented using Java 17, under our developing framework MORK.² In addition, all experiments have been executed on the same virtual machine, provided with 32 cores, 16 GB of RAM, and Ubuntu Server 20.04 as the operating system. We have also used Gurobi 9.1.2 as an exact solver in the replication of the state-of-the-art algorithms.

In the case of the DRFLP, state-of-the-art algorithms have been tested on a set of 38 synthetic instances. These instances have been proposed in several papers related to facility layout problems (Amaral, 2006, 2019; Anjos and Vannelli, 2008; Secchin and Amaral, 2019; Simmons, 1969). Additionally, we have extended the set of instances by introducing 15 new larger instances also used in the FLP problem family (Maadi et al., 2017).

According to the topology and structure of the DRFLP, instances can be understood as a weighted and undirected graph. The vertices correspond to facilities, which have a given size corresponding to the width of the facility. The weight of the edges will correspond to the pairwise material handling cost between facilities. If this value is zero, no edge will connect two vertices. Therefore, we apply some common graph metrics for both instance sets. Particularly, we show in Table 2 the average number of vertices (Avg. Vert.), the average number of edges (Avg. Edges), the average facility width (Avg. Width), the average edge weight (Avg. Weight), and density. The selected metrics are able to describe the main features of the instances, with no redundant information. Both instance sets and the analysis can be publicly found as artifacts as specified in Section 6.

Given the complete set of 53 instances, the proposed benchmark instance selection process was executed in order to determine the benchmark set. First, the general properties of the procedure are configured. Specifically, the number of instances of the benchmark set is established as 15% of the total sample size, that is, 8 instances. According to the number of components of the PCA, we select the minimum number of components whose explained variance is greater than 90%. Regarding the k-means algorithm configuration, k-means++ is chosen for the selection of the initial values of the algorithm (further details of the k-means++ initialization algorithm can be found in Arthur and Vassilvitskii, 2006). Besides, the maximum number of iterations is set to 1000 and the maximum number of clusters is set to 15. Once the input parameters for the process have been specified, the four steps detailed in Section 2.1 are performed.

The benchmark instance selection process starts with the structural features characterization. Considering the instances as graphs, we propose the following structural properties: number of vertices, number of edges, maximum/minimum/average/standard deviation of the degree of the vertices, maximum/minimum/average/standard deviation of the size of the vertices, maximum/minimum/average/standard deviation of the weight of the edges, and density of the graph. Hence, these 15 features were obtained for each one of the 53 instances.

Finally, the benchmark instance selection process concludes by selecting the most representative instances of each cluster, that is, the closest instances to the centroid of each cluster, following a descending order. The number of instances selected from each cluster depends on the cardinality of it. Specifically, the selected 8 instances are: 40-02, 40-04, 40-06, A60_03, A70_03, sko56_04, 14a, and Am13a. Moreover, in Table 3, we report the same graph metrics as in Table 2 for the selected instances and the cluster they belong to. As it can be easily observed, there is a clear distinction between instances of each cluster. For instance, selected instances from clusters 0 and 1, which share the same number of vertices, differ in other features like number of edges and average weight. Similarly, instances from cluster 3 present an average width very different than instances from cluster 2, which are closer in terms of number of edges. A similar analysis can be performed for all the different combinations of clusters.

Following the guidelines of the proposed methodology, the parameters of the algorithm were adjusted with an automated method. In particular, the algorithm has been tuned using irace (López-Ibáñez et al., 2016).

Table 4 summarizes the parameters given to irace. In particular, “constructive” decides which constructive implementation to use in the initial constructive method (step 1 in Algorithm 1). Their possible values are: “random,” corresponding to the random generator of solutions denoted as $C0$ in Section 4.2, “graspgr” and “grasprg,” which correspond to the greedy-random and random-greedy methods $C1$ and $C1'$ (Algorithm 2); and “tetris,” which corresponds to $C2$ (Algorithm 3). Parameter “alpha1” represents the $α1$ parameter for the initial constructive method. The reconstruction method (step 5 in Algorithm 1) is firstly configured by using the “reconstructive” parameter, which specifies the reconstructive implementation to use. Notice that it receives the same set of values as the initial constructive method. Parameter “alpha2” defines the $α2$ parameter for the reconstruction phase and “destratio” defines the percentage of the solution to destroy during the destruction phase, denoted as $β$ in Algorithm 1. Notice that there is no parameter for the destructive method (step 4 in Algorithm 1) since $D2$ destruction is required in case any variant of $C2$ is selected, since they both are based on combinations of partial solutions (see Section 4.4). On the contrary, $D1$ is selected if any variant of $C1$ is chosen by irace. The following parameters correspond to the multistart configuration. Parameter “stop” defines the strategy used to determine the threshold to stop (step 3 in Algorithm 4): by fraction or by constant values. Hence, depending on the selected value, a “fractionv”fraction for a given instance, or “constantv” constant value must be selected. The “increment” parameter defines how the fake facilities are incremented (step 7 in Algorithm 4): linearly with “linearinc” value or following the Fibonacci function with the “fiboinc” value. Parameter “linearratio” refines the value of the linear increment. Finally, since irace is able to choose the number of iterations for both the multistart and the iterated greedy components, irace will try to maximize both as they will never worsen the final result, getting an increment of the total computing time. Instead of providing a maximum number of iterations for both components to irace, the parameter “iterationsratio” is defined as the proportion of iterations spent in the outer for-loop of the MS-IG (⁠$tmax$ in Algorithm 4) compared with the inner for loop of the IG method (⁠$imax$ in Algorithm 1) such as outerIterations * IGIterations$=$ 1000000. This way, irace is able to balance the time spent diversifying with different solutions or intensifying working on the same solution.

This configuration is described in the “parameters.txt” file inside the irace folder in the repository specified in Section 6.

Using a maximum budget of 1000 experiments, and limiting the total execution time of any algorithm to 60 seconds, the best algorithm configurations found by irace are summarized in Table 5.

Elite configurations have several properties in common: first, the constructive chosen for all configurations is $C1$⁠. Both variants (greedy random and random greedy) have high values for “alpha1,” which means that the constructive method favors diverse solutions, accepting suboptimal moves during the construction phase. The best reconstruction method is $C2$⁠, and during reconstruction, the “alpha2” parameter is more restrictive. The preferred stop or threshold method is using a fraction of the instance size, and best results are achieved by incrementing the number of fake facilities linearly each iteration. Finally, the “iterationsratio” is balanced, or in other words, the number of iterations will be distributed uniformly between the iterated greedy and the multistart loops.

The best algorithm configuration found, represented as $#1$⁠, is selected for the comparison with the state of the art.

Once the best algorithm configuration has been determined by irace, its performance is assessed over the complete set of 53 problem instances.

Firstly, our multistart iterated greedy (MS-IG) proposal is compared with the four heuristics proposed in Amaral (2020) (H1, H2, H3, and H4) under the same set of 38 instances tackled in that paper. Table 6 illustrates this comparison, showing the averaged objective function value (Avg. o.f.), execution time in seconds, percentage deviation to the best value found (% Dev.), and number of best values found (# Best) for each algorithm, highlighting in bold font the best results. As seen, our MS-IG proposal obtains all the best values (38 of 38), spending short computation times. Detailed tables with objective function values and execution times per instance are available in the Appendix in Tables 8 and 9, respectively. Notice that a fair comparison of execution times cannot be made, since the previous work does not report the characteristics of the computer where the experiments were run.

Due to the fact that the previous work has not made any executable, artifact, or source code available, we have implemented their proposal following the descriptions in Amaral (2020). To ease future comparisons, we have made it available on the same code repository as our proposal, including an explanation of the correspondence between each coded instruction with Amaral's algorithm.

In this small replicability study, we have also assessed the validity of our implementation of the state-of-the-art method in order to verify that the comparison is fair. In this way, the results reported in Amaral (2020) are compared to our implementation of the previous methods using the original set of instances, obtaining a maximum, minimum, and average deviations over the results of $2.41%$⁠, $-0.02%$⁠, and $0.32%$⁠, respectively. Despite these differences, which could be attributed to fine-grain implementation details, we believe that our implementation is correct. Complete tables comparing both the objective function values and execution times between each heuristic pair for each instance can be found in the Appendix, in Tables 10 and 11, respectively. Notice that our implementations of Amaral's methods are labeled as H1', H2', H3', and H4'.

Note that, probably due to hardware differences and programming language used in the previous work computer and ours, the execution times are different. In particular, as seen in Table 11, our implementation of the state-of-the-art heuristics is up to 10 times faster.

Once the accuracy of the previous heuristics implementations has been established, a new experiment is executed with the set of 18 larger instances. The previous heuristics are limited to one-hour execution time, while the MS-IG is limited to a maximum time of 600 seconds. Both algorithms are executed in the same machine, using the same runtime libraries. A summary of the results is presented in Table 7, where the MS-IG proposal obtains all the best values, with an average execution time one order of magnitude lower. Notice that the previous heuristics reach the time limit in all runs. As with the previous experiment, detailed tables comparing both the objective function values and execution times can be found in the Appendix, in Tables 12 and 13, respectively. Hence, the MS-IG proposal is able to obtain the best results in all 53 problem instances, spending short computation times in relation to the state of the art.

As it can be seen Figure 8, our MS-IG approach reaches higher probability values than the state-of-the-art methods. In addition, the credible intervals of all the state-of-the-art variants overlap, which means that their behavior is comparable. On the other hand, the MS-IG approach reaches the highest probability of obtaining the best solution on average, 0.457, in relation to the other methods (⁠$<$ 0.15 on average). Moreover, the overall performance of this algorithm is clearly distinguished, since the credible intervals are not overlapped. Hence, the superiority of the proposed methods is not only supported by the results, but also for this Bayesian analysis.

All artifacts related to this article are publicly available at GitHub³ and archived in Zenodo,⁴ which in turn is indexed in OpenDOAR in order to ensure compliance with the FAIR principles (Wilkinson et al., 2016).

Following the recommendations from our methodology (see Section 2.3), the available artifacts are as follows:

• All instances used, both from the small and the new proposed set.

• Excel files with both the raw results data and the processed and summarized versions as used in this work.

• CSV and Excel files for the Bayesian analysis with the objective function results for all the instances and algorithms.

• Python Jupyter notebook used to load, analyze and generate the benchmark instance set. The notebook can be easily adapted to work with instances from other problems.

• Java source code for both our proposal and our implementation of the state-of-the-art algorithms.

• Documentation detailing how to execute and reproduce the results, including extending the algorithms or using a new set of instances.

• The Python Jupyter notebook used for benchmark instance selection.

• Small Docker scripts that can build, run, and deploy containers that easily run the proposal in any platform.

• The R script used to generate the Bayesian analysis shown in Figure 8.

• Generated diagrams and tables.

The proposed methodology as well as the provided artifacts ensure the reproducibility of our work on the DRFLP. Repeatability studies may be developed since all the executable artifacts and instances are publicly available. Reproducibility is possible since the source code for all the elements is available, allowing the researcher to modify any of the random factors. Replicability can be performed based not only on the provided source code, but also in the pseudocode descriptions of all the methods. Finally, generalizability is also possible since all the steps of the study in the DRFLP are detailed, and the source code is also provided. Therefore, any fixed factors could be modified to compare with the proposed work.

The aim of this article is two-fold. On the one hand, we propose a practical methodology to reproduce the work on an optimization problem solved by means of metaheuristics. On the other hand, we have improved the state of the art on the DRFLP applying the proposed methodology.

As seen in the article, not only have a set of precise steps been described on the methodology, that is, benchmark instances selection, parameter values tuning, and artifact generation, but we have also provided actual artifacts that can be used by researchers in any kind of reproducibility studies. In particular, a Python notebook implementing all the methods described in the benchmark instances selection is available. In addition, our optimization framework MORK is connected to irace to automatically tune algorithmic parameters. Also, the proposal is able to generate Docker containers with source code as well as including the same experimental environment of this work. Note that the applicability of these steps is not limited to using a specific software. Hence, a researcher may follow either the description in the article or the source code to develop their own implementation, tune the algorithmic parameters with a different external software, and produce their own corresponding artifacts.

Regarding the DRFLP, we propose the application of a multistart iterated greedy (MS-IG) algorithm, where four different constructive methods and two destructive methods have been proposed. We follow a model of the problem which introduces “fake facilities” as spaces between the actual facilities using two different strategies implemented in the MS-IG external loop. As stated before, we have followed our semiautomated methodology to select the instances and tune the algorithm in order to assure reproducibility studies.

In the experimental results, we have collected 38 instances from previous papers, and we also propose 15 new larger instances for this problem. In order to compare our MS-IG proposal in the larger instances, we have replicated the state-of-the-art method, obtaining similar results in the previous instances. This way, we have also shown the difficulties of these tasks when no artifacts are provided. As described in the experimental experience, our MS-IG proposal is able to obtain all 53 best solutions, spending a fraction of the execution time compared with the state of the art. Finally, a statistical analysis is conducted following a Bayesian approach in order to assess the performance of the different algorithms on the tackled set of instances.

Currently, we are studying different related lines regarding reproducibility, such as automatic characterization of instances. Furthermore, we pursue the creation of adimensional metrics to more precisely evaluate the performance of different implementations of algorithms over a problem. As in the case of this work, we will not only describe our findings, but we also will apply them on actual optimization problems.

This work has been partially supported by the Spanish Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033/FEDER, UE) under grant refs. PID2021-126605NB-I00 and PID2021-125709OA-C22; Spanish Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI/FEDER, UE) under grant ref. FPU19/04098; and Comunidad de Madrid y Fondos Estructurales de la Unión Europea with grant ref. P2018/TCS-4566.

In order to facilitate future comparisons, the complete set of results for every algorithm and instance combination is provided in this appendix.

Tables 8 and 9 show, respectively, the objective function values and total execution time in seconds obtained by the methods described in Amaral (2020) (denoted as H1 to H4) and our MS-IG proposal on the set of 38 instances from Amaral (2020). The best results are highlighted in bold font. As can be seen, the MS-IG proposal reaches all the best values, spending only a fraction of the total computing time spent by the previous methods.

Table 10 shows the results obtained by our implementation of the state-of-the-art heuristic, denoted as H1' to H4', in relation to the methods from Amaral (2020) for the set of 38 previous instances. The differences were calculated as a percentage deviation over the objective function score (columns %H1 to %H4). Total execution times in seconds are provided in Table 11.