The MOEADr Package: A Component-Based Framework for Multiobjective Evolutionary Algorithms Based on Decomposition

Multiobjective evolutionary algorithms based on decomposition (MOEA/D) represent a widely used class of population-based metaheuristics for the solution of multicriteria optimization problems. We introduce the MOEADr package, which oﬀers many of these variants as instantiations of a component-oriented framework. This approach contributes for easier reproducibility of existing MOEA/D variants from the literature, as well as for faster development and testing of new composite algorithms. The package oﬀers an standardized, modular implementation of MOEA/D based on this framework, which was designed aiming at providing researchers and practitioners with a standard way to discuss and express MOEA/D variants. In this paper we introduce the design principles behind the MOEADr package, as well as its current components. Three case studies are provided to illustrate the main aspects of the package.


Introduction
Multiobjective Optimization Problems (MOPs) [43] are problems in which multiple objective functions must be simultaneously optimized by the same set of parameters . MOPs are often characterized by a set of conflicting objective functions, which results in the existence of a set of optimal compromise solutions, instead of a single globally optimal one. In this way, multiobjective optimization algorithms are frequently characterized by their ability to find (representative samples of) this set of solutions with different compromises between the objective functions.
Multiobjective Evolutionary Algorithms based on Decomposition (MOEA/Ds), originally proposed by [65], represent a widely used class of population-based metaheuristics for solving MOPs [59]. MOEA/Ds approach the problem by decomposing it into a number of single-objective subproblems, which are then solved in parallel by a set of candidate solutions commonly referred to as a population.
Based on this general concept, a number of variations and improvements have been proposed over the past decade [59]. As is common in the general literature on evolutionary algorithms, many of these improvements are presented as monolithic entities, in which a fixed composition of operators and adjustments to the original algorithm is presented as a single, novel approach. Opposed to this monolithic approach, [5] have argued towards a more modular design of evolutionary algorithms in general, where an optimizer is seen as a composition of multiple, specialized components. This component-based approach allows researchers to clearly identify the contribution of each component, and facilitates the automated generation of new variants based on existing components, e.g., for the solution of specific problem classes. This approach also allows users to more easily test and implement new components, streamlining the development of new ideas and the reproducibility of results.
In this context, we have developed a component-oriented framework for MOEA/Ds, in which modules can be easily added, removed, modified or recombined by either users or automated testing and tuning programs. In particular, we defined a Variation Stack, which allows a very flexible way to model many different combinations and use cases of variation, local search, and solution repair operators.
This framework is implemented as the MOEADr package, which contains not only the original MOEA/D components, but several components found in more recent variations, all of which were recast as instantiations of the proposed component-based framework. The package is implemented in the R language and environment for statistical computing, and has been published on the Comprehensive R Archive Network (CRAN), at https://CRAN.R-project.org/ package=MOEADr. A development version is also available at the project repository, http://github.com/fcampelo/MOEADr. This paper describes the package and its underlying framework, and introduces the concepts necessary for a user to sucessfully apply MOEADr to their research or application problem. First, we provide a short primer on multiobjective optimization, as well as a short review of the MOEA/D (Section 2). We then review existing implementations of the MOEA/D, as well as existing implementations of MOP solvers in the R ecosystem, and locate our contribution in this context (Section 3).
Following that, we detail the framework of the MOEADr package, and present the MOEA/D components that it currently implements (Sections 4 and 5). This serves as a reference to the main features of the package, as well as a starting point for future contributions.
Finally, in Section 6 we describe the basic usage of the package with three case studies: A basic example of solving a benchmark MOP using a traditional MOEA/D algorithm; an example of using the MOEADr framework in conjunction with an automated algorithm assembling/tuning method to generate a new algorithmic configuration for a specific problem class; and an example of adding a custom component to the framework. Readers who are not concerned with the theoretical underpinnings of the package can skip straight to this section.
where n f is the number of objective functions, x ∈ R nv represents a candidate solution, n v is the number of decision variables, f (·) : R nv → R n f is a vector of objective functions, and Ω ⊂ R nv is the feasible decision space, such that where g i (·) : R nv → R, i = 1, . . . , n g and h j (·) : R nv → R, j = 1, . . . , n h represent the inequality and equality constraint functions, respectively. The image of the set Ω, denoted by f (Ω), defines the set of attainable objective values. Given two feasible solutions x i , x j ∈ Ω, we say that x i Pareto-dominates x j (written as f (x i ) ≺ f (x j ) or, equivalently, . . , n f } and f (x i ) = f (x j ). A solution x * ∈ Ω is considered Pareto-optimal if there exists no other solution x ∈ Ω such that f (x) ≺ f (x * ). The set of all Pareto-optimal solutions is known as the Pareto-optimal set, defined as PS = {x * ∈ Ω | x ∈ Ω : f (x) ≺ f (x * )}. The image of this set is referred to as the Pareto-optimal front, PF = {f (x * ) | x * ∈ PS}.
A widely used way to solve MOP using classical optimization methods is to represent the MOP as an arbitrary number of scalar optimization problems, which are built through aggregation functions such as the Weighted Sum (WS) or Weighted Tchebycheff (WT) [43] approaches. Each scalar optimization problem, generated by the aggregation function and a given weight vector, leads to problem in which the global optimum coincides with a particular Pareto-optimal solution of the original MOP. In this way, an estimate of the Pareto-optimal set can be obtained by solving a set of scalar aggregation functions. However, simply performing an independent optimization of these multiple scalar problems tends to result in difficulties for generating an adequate approximation of the Pareto set [43,17], particularly when a well-spread sample of the Pareto-optimal front is desired. [17] and [14] note that the original MOP described in (1) can be solved through a multiobjective evolutionary algorithm (MOEA), which is a populationbased approach that attempts to converge to an approximation of the Paretooptimal set in a single run. This feature enables a continuous exchange of information between the estimated solutions, which is useful to promote a proper approximation of the Pareto-optimal set. Among the algorithms commonly employed for the solution of continuous MOPs, we focus here on the class of MOEAs based on the explicit decomposition of the multiobjective optimization problem, which are briefly introduced below.

Multiobjective evolutionary algorithms based on decomposition
Multiobjective evolutionary algorithms based on decomposition (MOEA/Ds), originally proposed by [65], combine features of both MOEA approaches and classical scalarization approaches for tackling MOPs. In general, a MOEA/D decomposes a MOP, as defined in (1), into a finite number of scalar optimization subproblems. Each subproblem is defined by a weight vector, which is used in an scalar aggregation function to calculate the utility value of any solution for that particular subproblem. This set of subproblems is solved in parallel by iterating over a set of candidate solutions, commonly called the population. Aggregation functions and weight vectors are chosen so that (i) an optimal solution to a given subproblem is also Pareto-optimal for the original MOP; and (ii) the optimal solutions to the set of subproblems are well distributed in the space of objectives. It is thus expected that solving the subproblems provides a reasonable approximation of the Pareto-optimal front, regarding both convergence and diversity criteria.
Each subproblem has one incumbent solution that is directly associated to it, i.e., the population size is equal to the number of subproblems. At each iteration, one new candidate solution is generated for each subproblem, by the application of a sequence of variation operators to the existing population. This set of new solutions is compared against each incumbent solution based on their utility values on the respective subproblem. The best solution for each subproblem is maintained as its incumbent solution for the next iteration, following rules defined by a specific update strategy.
When generating new candidate solutions or comparing their performance against incumbent ones, the algorithm defines a neighborhood for each subproblem which limits the exchange of information between candidate solutions. This neighborhood provides a certain locality to the variation operators and update strategy, aiming to regulate the convergence speed and global exploration abilities of the algorithm.
Among the features that motivate the application of a MOEA/D to the solution of a MOP, [34] highlight a few that stand out in terms of their usefulness. First, it is generally simpler to handle objective value comparisons and, to a certain extent, diversity maintenance in a MOEA/D than in other MOEAs. This means that the MOEA/D is frequently able to return a set of reasonably well-spread (in the space of objectives) solutions, even when the number of subproblems is small or the number of objectives is high. Scaling techniques for attenuating the effects of objective functions with vastly different ranges are also easily incorporated into the MOEA/D structure, as are constraint-handling techniques.
Despite their specificities, we can characterize these different MOEA/D instantiations by defining the following component classes in the algorithm: • The decomposition strategy, which determines how the weight vectors are calculated and, consequently, how the MOP gets decomposed.
• The aggregation function, which uses the weight vectors to generate singleobjective subproblems to be solved.
• The objective scaling strategy, which defines how differences in the range of values of the objective functions are treated.
• The neighborhood assignment strategy, which determines the neighborhood relations among the subproblems. This strategy defines the locality of the exchange of information between candidate solutions when applying variation operators and updating strategies.
• A Variation Stack, composed of one or more variation operators, which generates new candidate solutions from the existing ones. Our general definition of a variation operator is a function that modifies a candidate solution in the space of decision variables, based on information about the problem structure and/or the current and past states of the population as a whole. Notice that this general definition also encompasses repair operators and local search operators as special cases of variation.
• The update strategy, which determines the candidate solutions to be maintained or discarded after each iteration, based on their utility values for specific subproblems and on their neighborhood relations.
• The constraint handling method, which defines how to treat points that violate problem constraints.
• The termination criteria, which determines when the algorithm stops the search and returns the set of solutions found.
Based on this decomposition of the algorithm into its individual components, it is possible to define a common framework from which specific MOEA/D variants can be instantiated. This framework is detailed in Section 4.

MOEA/D implementations
Several implementations of the original MOEA/D and some of its variants for continuous MOPs are available online, mainly in Java, C++ and MATLAB. 1 All implementations mentioned here have their source codes readily available for download from the links listed in the references. Table 3 summarizes these implementations. It is important to notice that there is no implementation native to R, and that while jMetal [46] and the MOEA Framework [27] represent object-oriented frameworks with some component-  The MOEADr package was motivated by a perceived need to facilitate not only the application of existing MOEA/D variants but also the development and investigation of new components, as well as the fast reproduction of newly published innovations based on a library of existing components, with minimal need for reimplementation.
Finally, it is worth mentioning that while there are no specific MOEA/D implementations native to R, a few packages implementing different, generalpurpose algorithms for multiobjective optimization are available. Table 3   Most of those packages offer a closed, individual algorithm, with the exception of Jakob Bossek's ecr package [8], which offers a modular approach for instantiating evolutionary algorithms for both single and multiobjective optimization. For the latter class, the algorithmic structure defined by this package lends itself easily for the implementation of dominance and indicator-based approaches, but not necessarily for decomposition-based algorithms such the MOEA/D. This is also the case of another recent initiative to define a common framework for defining multiobjective evolutionary algorithms [4]. To the authors' knowledge, there is no R package implementing specific MOEA/D algorithms, nor any component-based implementation that allows an easy instantiation of multiobjective evolutionary algorithms based on decomposition.

The MOEADr package
The MOEADr package is an R implementation of multiobjective evolutionary algorithms based on decomposition, including many MOEA/D variants found in the literature. Our goals in this package are three-fold: • To provide a component-wise perspective on the MOEA/D family of algorithms, where each different algorithm exists as a configuration of common components; • To be easily extensible, so that researchers and users can implement their own components, facilitating both applied uses and scientific inquiries in this field; • To include a representative section of the existing literature in MOEA/D variants; To achieve these goals, the implementation was guided by the following design decisions: 1. For the sake of uniformity in the implementation of each component, we define the main variables of the MOEA/D: the solution set, neighborhood set, subproblem weight set, subproblem utility value set and violation value set as matrices; 2. Algorithms in the MOEA/D family are broken down into common components, and each component is recast as an operator on the matrices defined above; 3. Each individual component avoids, to the maximum degree possible, to produce and rely on side effects beyond the explicit manipulation of the main variables above; 4. Each component is implemented as a separate R script file with a fixed naming scheme. The choice of components is specified as a parameter to the main function call; 5. In particular, we define a Variation Stack, which is a list of variation components to be used by the MOEA/D in order. The choice of variation operators shows a very large diversity in MOEA/D literature, and using a variation stack allows for a uniform description of the many existing configurations; The general framework of the MOEA/D, as implemented in our package, is presented in Algorithm 1. The user defines an instance of this framework by choosing a specific component for each of the roles described in section 2.1, and one or more variation components to compose the variation stack V.
In the following sections we detail each of these component classes, presenting a formal definition of their roles as well as relevant examples. The list of components described in this paper (and currently available in the MOEADr package) is summarized in Table 4.

MOEA/D components available in the MOEADr package
In this section we describe in detail each of the components included in the MOEADr package, version 2.1. We present a formal definition of their roles as Constraint handling Table 3: Components currently available in the MOEADr package.

7:
Copy incumbent solution set X (t) into X (t)

13:
Update run flag Stop criteria (Sec. 5.7) 14: well as relevant examples, both from the specific MOEA/D literature and the wider body of knowledge on multiobjective evolutionary algorithms in general. The complete list of roles and specific components is summarized in Table 4.
For each component, the different notations of existing works were recast to a standard mathematical notation, to prevent confusion and ambiguities that would inevitably arise if we tried to follow the many different nomenclatures from the literature. Whenever possible, our notation employs vector and matrix operations to describe the modules, in an attempt to highlight the mathematical structure of each component, as well as differences and similarities among variants.
Special care was taken to guarantee the modularity of the definitions provided in the following sections, so that each component is independent from design choices made for the others. This characteristic allows the free exchange of components while guaranteeing the correct flow of the MOEA/D, at the cost of some implementation overhead. This also simplifies the use of automated algorithm assembly and tuning methods, as well as efforts for replicating and testing algorithms from the literature.
To describe the components, let the following common variables be defined: n f is the number of objective functions that compose the MOP. N ∈ Z >0 is the number of subproblems, which is either a user-defined parameter, or calculated internally by the decomposition strategy.
denotes the set of incumbent solutions for the subproblems at iteration t. With some abuse of notation, X (t) will also denote a matrix with each row i defined by is the matrix of subproblem weights calculated by the decomposition strategy (Section 5.1), and B ∈ Z N ×T >0 is the matrix of neighborhood relations calculated by the neighborhood assignment strategy (Section 5.3), with T denoting the neighborhood size.

Decomposition strategies
A decomposition strategy is responsible for generating the weight vectors that characterize the set of scalar subproblems in a MOEA/D. The number of subproblems can be either explicitly defined, or calculated indirectly by the decomposition strategy based on user-defined parameters.
Methods for generating aggregation weight vectors in the structure of MOEA/D are usually borrowed from the theory of design and modeling in experiments with mixtures [9]. Among the main designs in this area, the simplex-lattice design [55] has been the most commonly adopted in the context of decompositionbased MOEAs. This approach, as well as a multiple-layer variant are presented below. A third approach, based on uniform designs, is also described.
The following definitions hold for all strategies discussed in this section: let be a matrix of non-negative values. Also, let each row of Λ be a vector summing to unity, λ i 1 = 1. Each row λ i can be interpreted as the weight vector defining the i-th subproblem, and the matrix Λ is referred to as the weight matrix. Since, given a scalar aggregation function (Sec. 5.2), each weight vector defines a unique subproblem, λ i will also be used to refer to the i-th subproblem in the following sections.
In the MOEADr package, a list of available decomposition strategies can be generated using the function get decomposition methods().

Simplex-lattice design (SLD)
In the simplex-lattice design [9,65] the user provides a parameter h ∈ Z >0 that defines both the number of subproblems and the values of the weights. In this method, elements of the weight matrix can only assume h + 1 distinct values The simplex-lattice design generates N weight vectors by using all combinations (with repetition) of n f values taken from the set of values defined in (3), which results in subproblems defined when this strategy is used. In this way, a general (h, n f )simplex-lattice can be used to represent N weight vectors on the objective domain. For instance, for an arbitrary three-objective problem (n f = 3), a value h = 18 yields N = 20 2 = 190 distinct weight vectors.

Multiple-layer simplex-lattice design (MSLD)
To obtain a reasonable distribution of weight vectors within the n f -dimensional simplex using the simplex-lattice design it is necessary that h ≥ n f [20]. While this condition is generally harmless for MOPs with few objectives, a large number of weight vectors is generated for high-dimensional objective domains, even in the limit condition h = n f . On the other hand, making h < n f results in weight vectors sparsely distributed only at the boundary of the simplex, which jeopardizes the exploration abilities of the algorithm.
To address this issue, a multiple-layer approach can be devised. This strategy generates k subsets of weight vectors by generalizing the user parameter h to a vector of positive integers h = (h 1 , . . . , h K ) , h k ∈ Z >0 ∀k. Each constant h k is used to calculate N k weight vectors according to the SLD method and, within each subset, the vectors are scaled down by a factor τ k ∈ (0, 1] using a simple coordinate transformation with λ ki being the i-th vector in the k-th subset. Each subset must be associated with a unique user-defined value of τ k . 4 The total number of subproblems defined in this method is given by N = K k=1 N k , and the final weight matrix Λ is composed by the scaled weight vectors from all layers. This method ultimately amounts to generating weight vectors in layers, with each layer corresponding to a different-sized simplex in the space of objectives. This approach trades the loss of exploration ability by the increased number of user-defined parameters (it requires the definition of 2K + 1 parameter values). Common sense seems to suggest τ k = k/K as a reasonable heuristic for defining the scaling factors, but there is currently neither empirical nor theoretical support for this choice.
Finally, it must be remarked that this approach to generating the weight vectors generalizes the method of [35], which was defined specifically for K = 2.

Uniform design (UD)
The uniform design method represents an alternative approach to generate the weight vectors. Its use was proposed by [58], with the stated objectives of improving the distribution of the weight vectors and providing greater control over the number of subproblems.
The UD method for calculating the weight vectors can be described as follows. First, let H N be defined as the set where gcd (·) returns the greatest common divisor of two integers, and N ∈ Z >0 is a user-defined value that determines the number of subproblems. Let h = h 1 , . . . , h n f −1 be a vector composed of n f −1 mutually exclusive elements of H N . Any such vector can define a matrix U N (h) ∈ (Z >0 ) N ×(n f −1) with elements calculated as u ij = (ih j ) mod N . Denoting the set of all possible h vectors defined from H N as P (H N ), the next step is to identify the vector that results in the U N (h) matrix with the lowest CD 2 discrepancy where CD 2 (·) is the centered L 2 -discrepancy of a matrix [21,58], calculated as If we define U N = U N (h) − 0.5 /N , then the elements of the weight matrix Λ are returned by the transformation with the guarantee that Λ conforms with the properties stated at the beginning of section 5.1.

Scalar aggregation functions
The scalar aggregation function is used to calculate the utility value of candidate solutions for each subproblem. This is done by specific functions of the objective function values f (x) and the weight matrix Λ.
In the MOEADr package, a list of available aggregation functions can be obtained using function get scalarization methods().

Weighted sum (WS)
This technique performs a convex combination of the objective values, resulting in scalar problems of the form where z = (z 1 , . . . , z n f ) is a reference vector with the property that ∀j z j ≤ min{f j (x) | x ∈ PS}. A diverse set of Pareto-optimal solutions can be achieved by using different weight vectors λ i . However, due to the convex nature of this operation, only convex sections of Pareto fronts can be approximated using this strategy [43].
It is important to highlight here that obtaining a good estimate of z is by itself a computationally-intensive effort. In many cases, it is common practice to iteratively estimate this reference point from the set of all points visited up to a given iteration, In these cases the elements of the estimated reference vector,ẑ (t) , are defined aŝ This estimated vector is updated at each iteration and is frequently employed instead of a fixed z in the scalar aggregation functions used with the MOEA/D.

Weighted Tchebycheff (WT)
In this approach, the scalar optimization problems are defined as where denotes the Hadamard product, and · ∞ is the Tchebycheff norm. Unlike the weighted sum approach, the weighted Tchebycheff approach is not sensitive to the convexity of the Pareto front [43]. However, this approach offers poor diversity control, as the solution of (11) for equally-spaced weight vectors do not necessarily translate to a well-spread approximation of the Pareto front [51,61].

Penalty-based boundary intersection (PBI)
The PBI aggregation function [65] is an extension of an earlier approach known as normal boundary intersection [15]. The scalar optimization problems defined by the PBI strategy are given as where · 2 denotes the Euclidean norm, and θ pbi ∈ R ≥0 is a user-defined penalty parameter. Notice that d i,1 is related to the convergence of f (x i ) towards the Pareto-optimal front, whereas the minimization of d i,2 provides a way to control solution diversity. This aggregation function enables the definition of a trade-off between convergence and diversity (in the space of objectives) by an a priori adjustment of θ pbi , which directly influences the performance of the MOEA/D.

Adjusted weighted Tchebycheff (AWT)
An alternative that attempts to address the poor diversity control of the weighted Tchebycheff approach is the adjusted (or transformed) Tchebycheff scalarization function [51,61]. This method defines the scalar subproblems as: with the elements of the vector ρ i defined as where is a small constant added to prevent divisions by zero. 5 Notice that the only difference between (11) and (13) is the substitution of the weight vector λ i by its respective "normalized inverse" ρ i . Besides minimizing the distance between f (x) and z, this transformation promotes the distance minimization between an objective vector f (x) and its corresponding vector ρ i [51]. In this aspect, the adjusted weighted Tchebycheff strategy presents an idea similar to that of the PBI, but without any additional parameters.

Inverted PBI (iPBI)
While the PBI approach defines its search based on an ideal reference point z which represents an (estimated) ideal solution, the inverted PBI function uses a nadir reference point z, defined as ∀k z k ≥ max{f k (x) | x ∈ PS} [53,54]. 6 This method works by defining the scalar minimization 7 problems with where, similarly to the PBI approach, θ ipbi ∈ R ≥0 is a user-defined parameter that controls the balance between d i,1 (convergence) and d i,2 (diversity).

Scaling of the objective domain
Since the range of objective functions in MOPs can present arbitrarily large differences, an appropriate scaling of the objective domain is sometimes used to improve the performance of the MOEA/D. When performed, this scaling happens prior to the calculation of the scalar aggregation function, so that scaled function valuesf i (·) replace f i (·) in the calculations. Let k ) ∀k = 1, . . . , N be estimates of the ideal and nadir objective vectors at a given iteration t. A straightforward way to standardize the objective domain [43,65] is to replace each function value f i (x) bȳ which guarantees thatf i (x) ∈ [0, 1] ∀i. 8

Neighborhood assignment function
The neighborhood assignment function generates a matrix B defining the neighborhood relationships between subproblems. In general, neighborhood relations among the subproblems are used for defining restrictions to the exchange of information among candidate solutions when applying the variation operators, as well as for regulating the replacement of points at the end of every iteration [29]. Neighborhood relations can be defined either in the space of decision variables, Ω, or in the space of objectives f (Ω). The more usual case is the definition of neighborhoods based on the distances between weight vectors in the space of objectives [65,34]. Let M ∈ R N ×N ≥0 be the symmetric matrix defined by taking all pairwise Euclidean distances between weight vectors, i.e., with elements given by Let m i denote the i-th row of M. For the i-th subproblem, its neighborhood vector b i ∈ Z T >0 consists of the indices of the T smallest elements of m i , with T ∈ Z >0 a user-defined parameter. 9 Since weight vectors are usually kept constant throughout the optimization process (with some exceptions, see [51]), the neighborhoods have to be determined only once in this approach, and remain fixed throughout the execution of the algorithm. An alternative approach uses Euclidean distances between incumbent solutions in the space of decision variables [12] to define the neighborhood relations among subproblems. In this case the distance matrix M is defined by the distances between all vector pairs x with x (t) i ∈ Ω being the incumbent solution to the i-th subproblem at iteration t. A neighborhood vector b i is then composed by the indices of the T subproblems whose incumbent solutions are closest, in the space of decision variables, to the one associated with λ i . As the incumbent solutions change across iterations, neighborhood relations must be updated whenever a new incumbent solution is determined for any subproblem. While this results in increased computational cost (N (N − 1)/2 distance calculations per iteration), this neighborhood definition may contribute to the algorithm performance in problems for which solution similarity in Ω does not correspond to small distances in f (Ω) [12].

Variation stack
As in any evolutionary algorithm, variation operators in the MOEA/D generate new candidate solutions based on information about points visited by the algorithm up to a given iteration. The standard MOEA/D [65] employs two variation operators, namely Simulated Binary recombination (SBX) [19] followed by Polynomial Mutation [18], a combination still widely used in the literature [2,35]. Another very successful MOEA/D version known as MOEA/D-DE [34,58] employs Differential Mutation and Binomial Recombination [57], followed by Polynomial Mutation.
While these two sets of variation operators are arguably the most widespread in the literature, any combination of variation operators can, at least in principle, be used to drive the search mechanism of the MOEA/D, as long as they are sequentially compatible. This includes the possibility of employing multiple combinations of recombination and mutation variants -applied either sequentially or probabilistically -as well as the use of local search operators. Therefore, Algorithm 1 employs a user-defined set of variation operators which are applied sequentially to a population matrix X (t) , copied from the incumbent solution matrix X (t) prior to the application of the variation operators.
In the following descriptions all operators (with the exception of the repair operators described later) are defined assuming that the decision variables are contained in the interval [0, 1], which greatly simplifies many calculations.
Whenever a MOP is defined with decision variables bound to different limits, it is assumed that the variables are properly scaled prior to the optimization.
Also included in this discussion of variation operators are Local Search Operators, which execute a limited search of the solution space focused on the region close to a particular solution. Because they act as an operator that directly modifies the solution set, in this framework we group them together with other variation operators.
The variation operators available in the MOEADr package can be listed using functions get variation operators() and get localsearch methods().

Definition of sampling probabilities for variation
Prior to the application of the variation stack, the neighborhood effects (Sec. 5.3) must be determined. For each subproblem, a vector p i = {p i,1 , p i,2 , . . . , p i,N } is defined as containing the probabilities of sampling the candidate solutions associated with each subproblem when variation operators are applied for the i-th subproblem.
In a general form, sampling probabilities are defined by the neighborhood vectors b i and by a user-defined parameter δ p ∈ [0, 1], which represents the total probability of sampling from the specific neighborhood b i . Probabilities p i,j are defined as The standard MOEA/D [65] samples exclusively from the neighborhood of each subproblem (δ p = 1 and, consequently, p i,j = 0 ∀j / ∈ b i ). This setting intensifies local exploration, at the cost of a loss of solution diversity, which may compromise the effective exploration of the design space in later iterations [34]. Less restrictive approaches have been used in the literature [34,12], defining δ p < 1, albeit usually at relatively high values. In these approaches, while each subproblem maintains a strong bias towards using information from subproblems indexed by b i , access to the other candidate solutions remains possible, enabling the generation of a more diverse set of points by the variation operators.

SBX recombination
Let η X ∈ R >0 be a user-defined parameter, and u i ∈ [0, 1] nv be a vector of uniformly distributed random values. Also, let β i ∈ R nv be a vector with elements defined as Let bi ∈ R nv be two vectors sampled from X (t) according to the sampling probabilities defined for the i-th subproblem (18). For each subproblem, SBX recombination [19,65] produces a new candidate solution according to where u i ∈ [0, 1] is a uniformly distributed random value, and p X ∈ [0, 1] is a user-defined parameter. After this operator is applied for all i ∈ {1, . . . , N }, each row of the population matrix X (t) is updated with its corresponding point x (t) i . 10

Polynomial mutation
Let η M ∈ R >0 and p m ∈ [0, 1] be user-defined parameters, and u i ∈ [0, 1] nv be a vector of uniformly distributed random values. For a given candidate solution x (t) i ∈ X (t) , let β i ∈ R nv be a vector defined by if u i,j ≤ 0.5, or otherwise, with η = η M + 1. Let also v i ∈ {0, 1} nv be a vector with elements sampled independently from a Bernoulli distribution with probability parameter p M . For each subproblem, the polynomial mutation [18,65] generates a candidate solution according to After this operator is applied for all i ∈ {1, . . . , N }, X (t) is updated with the new points x (t) i .

Differential mutation
Let φ ∈ R =0 be a user-defined parameter 11 , and x (t) bi ∈ R nv be two mutually exclusive vectors sampled from X (t) according to the sampling probabilities 10 Despite being much more compact than the usual definitions of SBX recombination [19,65], this one is equivalent provided that the stated scaling of the variables to the interval x i,j ∈ [0, 1], ∀i, j is maintained. 11 φ can be set either as a constant value, or defined randomly for each application of this operator. In the latter, it is common to independently sample φ ∈ (0, 1] for each operation, according to a uniform distribution [39]. defined in (18). For the i-th subproblem the differential mutation operator [50] generates a new candidate solution as where x (t) i,basis ∈ R nv is a basis vector, which can be generated in several ways. The most common strategy used with the MOEA/D [34] is to randomly sample a vector from X (t) (mutually exclusive with x (t) ai and x (t) bi ) according to the sampling probabilities defined in (18).
Two possible alternatives for the basis vector in the context of MOEA/D are also suggested here. For a given subproblem i, let the points x be ordered in decreasing order of utility (i.e., in increasing order of aggregation function value). One possibility is to use the mean point of the neighborhood while the second is to use a weighted global intermediate point, which is an adaptation of the basis vector commonly used in evolution strategies [1] x (t) with weights given by Further alternative definitions of the differential mutation operator in the MOEA/D context can be found in the literature [39]. Regardless of the specifics of each definition, this operator is known to be invariant to rotation [50,39]. After its application, X (t) is updated with the points x (t) i .

Binomial recombination
Let ρ ∈ [0, 1] be a user-defined parameter, k i ∈ {1, . . . , n v }, i = 1 . . . , N denote a set of randomly selected integers, and u i ∈ [0, 1] nv be a vector of uniformly distributed random values. Also, recall that the incumbent solutions at iteration t (unmodified by any variation operators) are stored in the matrix X (t) . This operator can then be expressed by the sequential application of and (29) x (t) As with the previous operators, after this one is applied to all subproblems, each row of X (t) is updated with its corresponding point x (t) i .

Repair operators
Repair operators in evolutionary algorithms for continuous optimization usually refer to strategies for ensuring that the box constraints, i.e., pairs of constraints defined by (x j,min − x i,j ≤ 0 ; x i,j − x j,max ≤ 0) ∀j, are satisfied. While the literature tends to place repair operators in a separate category from the variation ones, we argue that they actually belong to the same class of operations in the MOEA/D framework, as both are used to modify candidate solutions according to specific rules.
While there are a number of repair methods [28], we describe here only the simple truncation repair, which can be defined as where functions max(·, ·) and min(·, ·) return the largest and the smallest of their arguments, respectively. After this operation is performed, X (t) is updated accordingly.

Local search operators
Local search strategies are usually employed in the MOEA/D to accelerate convergence. As with the repair operators, these approaches are commonly classified as distinct from the variation operators, but we argue that they are indeed part of the same block, for the same reasons stated earlier.
Since these operators sometimes result in the need for additional evaluations of candidate solutions, as well as in loss of diversity, their application is usually regulated by a frequency parameter, expressed either as a period of application τ ls ∈ Z >0 (i.e., the operator is applied to each subproblem once every τ ls iterations); or as a probability of application γ ls ∈ [0, 1] (i.e., at every iteration the operator is applied with probability γ ls for a given subproblem) [17,14]. Whenever any of these conditions is true, local search is applied on the incumbent solution of a subproblem instead of all other variation operators. Notice that both τ ls = 1 and γ ls = 1 lead to the same behavior, that is, the MOEA/D ignoring all other variation operators and performing only local search at every iteration. Values that are sometimes offered for these parameters in the wider MOEA literature are τ ls = 5 [63] and γ ls ∈ [0.01, 0.05] [17,14], albeit without much theoretical or experimental justification.
As with all other variation methods, after this operation is performed the rows of X (t) are updated with the corresponding points x ∈ R nv be the three candidate solutions with the highest utility for the ith subproblem at iteration t, such The three-point quadratic approximation method [58] for performing local search in the MOEA/D can be expressed as where where ∈ R >0 is a small positive constant. 12 Differential vector-based local search (DVLS): Although this strategy was originally proposed for non-decomposition MOEAs [10], its adaptation to the MOEA/D framework is straightforward. Let ∈ X (t) be two candidate solutions, with (a i , b i ) denoting mutually exclusive indices sampled from the neighborhood b i . Two new candidate alternatives, x + i and x − i , are then generated according to (24), using x and φ = ±φ ls , with φ ls ∼ N (µ = 0.5, σ = 0.1) [10]. This local search operator then compares these two alternatives against the incumbent solution, and returns the point with the best performance for the subproblem, i.e., 12 [58] recommend = 10 −6 .

Update strategies
Update strategies in the MOEA/D play the same role as selection operators in general MOEAs, regulating the substitution of incumbent solutions by those resulting from the application of a sequence of variation operators. In the MOEA/D, this substitution is controlled by the replacement strategy [65,2,62], as well as by the neighborhood relations between subproblems. Three update methods are presented next.
In what follows, let X (t) denote the matrix containing the candidate solutions generated after the application of the variation stack, and X (t) be the matrix containing the incumbent solutions at iteration t. Also, let f agg (·) denote the aggregation function (see Sec. 5.2). For each subproblem i, the update strategy will generate a set of candidate solutions C i ( generally composed of the incumbent solution x The list of update strategies available in the MOEADr package can be generated using get update methods().

Standard neighborhood replacement
Let the candidate set for the i-th subproblem be defined as where b i is the neighborhood defined in Section 5.3. This replacement strategy [65] updates the population according to (37) That is, the best solution to the i-th subproblem considering the incumbent solution and the candidate solutions indexed in b i .

Restricted neighborhood replacement
In the standard replacement, a single candidate solution can replace up to T (the neighborhood size) incumbent solutions at any given iteration, which can lead to diversity loss and premature stagnation. To overcome this issue, [34] proposed a limit to the number of copies that any single point x can pass on to X (t+1) , defined by a user-defined parameter n r ∈ Z >0 | n r ≤ N . This approach generalizes the standard neighborhood replacement (which happens if n r = T ), and provides a relatively simple way to limit the propagation speed of candidate solutions in the population. The definition of good values for n r , however, may require some tuning: both n r = 2 [34] and n r = 0.01N [66] are recommended, without much theoretical or empirical support.

Best subproblem replacement
In the previous strategies, a candidate solution x (t) i is considered for replacing the incumbent solutions of its neighboring subproblems, i.e., it can replace incumbent solutions x (t) k | k ∈ b i . However, it may happen that the best available candidate solution for a given subproblem i is not among those in the neighborhood b i . To address this issue, [62] suggested a three-step replacement strategy. The first step is to identify the subproblem k i for which each new candidate solution x (t) i is most effective, i.e., The second step is to generate replacement neighborhoods b r ki , which are composed of the T r nearest neighbors to the k i -th subproblem (calculated as in Sec. 5.3), where T r ∈ Z >0 | T r ≤ N is a user-defined parameter. Finally, for each subproblem i, the original neighborhood b i is replaced by b r ki , and then the Restricted Neighborhood Replacement is used for updating the population.

Constraint handling approaches
Some possible approaches to deal with constraints in the MOP formulation are discussed below. While most available implementations of the MOEA/D (see Section 3) provide versions that are only capable of dealing with box constraints (using repair operators such as truncation), dealing with general constraints in the context of the MOEA/D is relatively straightforward, as discussed below.
The constraint handling techniques available in the MOEADr package can be consulted using get constraint methods().

Penalty functions
The most common approach in the MOEA community to handle constraints is to use penalties [43,13]. The idea of penalty functions is to transform a constrained optimization problem to an unconstrained one, by adding a certain value to the aggregation function value of a given point based on the magnitude of constraint violation that it presents.
Assuming that the magnitude of the constraint violation is derived from the sum of violations of all inequality and equality constraints, its value can be defined as 13 where ∈ R ≥0 is a small threshold representing a tolerance for the equality constraints. The new (penalized) aggregation function to be optimized is then obtained as f agg in which β v ∈ R >0 is a user-defined penalization constant. For each subproblem, the penalized values are then used to select which solution will become (or remain) the incumbent one in the next population, X (t+1) , according to the update strategies (Sec. 5.5).

Violation-based ranking (VBR)
This constraint handling technique generalizes a few methods available in the literature, such as Tournament Selection (TS) [16], Stochastic Ranking (SR) [52], and Violation Threshold (VT) [3]. Like the Penalty Functions approach, Violation-based Ranking uses the total magnitude of constraint violations v(x) (39) to penalize unfeasible solutions, but this penalization takes the form of ranking functions that employ different quantities depending on the feasibility (or not) of the solutions being compared.
In its general form, this constraint handling method uses the following criteria when ranking the solutions: • If a solution is feasible, use its f agg (x|λ k , z) value; • If a solution is unfeasible, check a criterion c (x): -If c (x) = T RU E, use its f agg (x|λ k , z) value; -If c (x) = F ALSE, use its v (x) value; • Assuming that n 1 ≤ N points are to be ranked using their f agg (x|λ k , z) values, VBR ranks those solutions first (rank values from 1 to n 1 ), and then attributes ranks for the remaining solutions using their v (x) values, from n 1 + 1 to N . Rank values are then used instead of f agg (x|λ k , z) whenever solutions must be compared.
As mentioned earlier, this criterion generalizes the commonly-used constraint handling techniques Tournament Selection (TS), Stochastic Ranking (SR), and Violation Threshold (VT). These specific methods can be instantiated from the rules defined in Table 5.6.2, by choosing an appropriate function c (x), see Table  5.6.2. For SR, p f ∈ [0, 1] represents a user-defined parameter, and u ∈ [0, 1] denotes a uniformly distributed random value. For VT, the adaptive threshold v is calculated independently for each subproblem, as
An important point to highlight here is that the VBR method does not necessarily guarantee that feasible solutions will always be maintained in the population. This means that a good feasible solution obtained previously for a scalar subproblem may be eliminated later in the search. To prevent this loss of feasible solutions, an elitist archiving strategy can be adopted in order to preserve the feasible solution with the best aggregation value found for each sub-problem throughout the iterations. In this case, this elitist archive should be considered as the output population of the algorithm [64].
Finally, it is worth mentioning that this generalization of known constraint handling techniques may serve as a template for further developments of better methods for dealing with constraints in population-based algorithms, both for single and multi-objective cases. This possibility, however, falls outside the scope of the current paper, and will not be explored further.

Termination criteria
Termination criteria determine when the algorithm should stop running and return the solutions found. Common criteria in the evolutionary computation literature tend to fall into two broad categories, namely: (i) exhaustion of the algorithmic budget, in terms of execution time, number of points visited, or number of iterations performed; and (ii) attainment of a specific threshold value for some criterion, such as lack of improvement, loss of population diversity, convergence to specified quality values, or some other population statistic. 15 Time-based criteria stop the process after a user-defined amount of time has passed since the program started running. These criteria can be useful when comparing computationally intensive tasks such as the construction of neighborhood or subproblem weight matrices, but it can be sensitive to background tasks if wall clock is used. In general, process time tends to be preferable, although it tends to add some complexity to the implementation.
Number of iterations criteria stop the search immediately after the iteration counter becomes higher than a given threshold. This tends to be useful when comparing minor variations of a given algorithm, as long as there is no significant differences in the computational cost per iteration.
Finally, number of evaluations criteria terminate the search when the total number of points visited by the algorithm reaches a predefined threshold. This criterion is usually preferred in situations when the larger part of the total computational cost is incurred by utility function evaluations, in which case it is a good practice to keep this value constant when comparing very different algorithms. Notice that because the termination check occurs once per iteration, the total number of evaluations may show some fluctuations from the userdefined value.
The list of available stop criteria in the MOEADr package can be consulted using function get stop criteria().

Usage examples
In this section we present three case studies which illustrate common usage situations for the MOEADr package. The first shows how to solve a simple MOP using the package. The second illustrates the application of the component oriented framework with an automatic algorithmic assembling and tuning approach. Finally, the third example demonstrates how to use a user-defined component with the package.
These three case studies are also available as vignettes in the MOEADr documentation.

Solving a simple multiobjective optimization problem using MOEADr
In this case study, we use MOEADr to optimize a simple 2-objective problem composed of the (shifted) Sphere and Rastrigin functions in the domain Ω = [−1, 1] nv , defined as The R implementation of the problem can be defined as follows. Note that the MOEADr package requires the multiobjective problem to be defined as a func-tion that can receive a population matrix X, and return a matrix of objective function values for each point. The moead() function requires a problem definition, discussed above, an algorithm configuration, logging parameters, and a seed. In this example we used an algorithm preset. The preset moead() function can output a number of different presets based on combinations found on the literature. These presets can also be modified (either partially or as a whole), as will be shown in another case study. preset moead("original") returns the original MOEA/D configuration, as proposed by [65]. Running preset moead() without an argument outputs a list of available presets.
The moead() function returns a list object of class moead, containing the final solution set, objective values for each solution, and other information about the optimization process. The MOEADr package uses S3 to implement versions of plot(), print() and summary() for this object.
plot() will show the estimated Pareto front for the objectives, as in Figure  1. When the number of objectives is greater than 2, a parallel coordinates plot is also produced (see Figure 2). summary() displays information about the number of non-dominated and feasible solution points, the estimated ideal and nadir values, and (optionally) the inverted generational distance (IGD) and hypervolume indicators [68]

A more complex example
Package smoof [7] provides the implementation of a large number of single and multiobjective test functions. MOEADr provides a wrapper function make vectorized smoof() to easily convert smoof functions to the format required by the moead() function, as illustrated below for a 5-objective standard MOP benchmark problem. The code below shows an example of how to modify an algorithm preset. Because of the higher number of objectives, we want to reduce the value of the parameter H in the SLD decomposition component (see Section 5.1) used by the preset from 100 to 8: As before, we see an overview of the results with summary() and plot(). Notice that the output of plot() is different when the number of objectives in a problem is greater than 2, as shown in Figure 2.

Fine tuning algorithm configurations using MOEADr and irace
For this example, we employ the Iterated Racing procedure (available in the irace package) [41], to automatically assemble and fine-tune a MOEA/D configuration based on the components available in the MOEADr package.
Ten unconstrained test problems from the CEC2009 competition 16 are used, with dimensions ranging from 10 to 30. Dimension 20 was reserved for testing, while all others were used for the training effort. To quantify the quality of the set of solutions returned by a candidate configuration we use the Inverted Generational Distance (IGD) [68] indicator. The number of subproblems was fixed as 100 for n f = 2 and 150 for n f = 3.
We define a tuning budget of 1,000 runs for this example. The algorithm was set with the following fixed parameters: Uniform decomposition, AWT agreggation function, restricted update, simple scaling, max. evaluations = 50, 000 as a stop criterion. The variation stack was partially fixed with the sequence of operators differential mutation (with φ ∼ U (0, 1)), binomial recombination, polynomial mutation (with p M = 1/n v ) and simple truncation repair.
Seven tunable parameters were set as follows: • Neighborhood type ∈ by λ i ; by x (t) i ; • Neighborhood size T ∈ [10, 40]; • Probability of sampling from neighborhood δ p ∈ (0.1, 1.0); • For the restricted update: n r ∈ [1, 10]; • Differential mutation basis ∈ {rand, mean, wgi}; • Binomial recombination ρ ∈ (0, 1); • Polynomial mutation η M ∈ (1, 100); The code of this case study requires a large amount of boilerplate for the irace package [41] and, for the sake of brevity, it is not included in the text. Please refer to the supplementary materials for the full replication script. 17 Figure 3 shows the IGD values achieved by the four final configurations over the test problems. Based on these final configurations, we assemble the consensus MOEA/D configuration, which is described in Table 6 Table 6: Final MOEA/D configuration returned by the iterated racing procedure. "Value" shows the best configuration, and "Consensus" represents the agreement between the four final configurations.
literature. The example used in this paper is intended only as a proof-of-concept, but we highly recommend that a similar approach is used when developing new components, in order to observe not only the individual performance of the novel component, but also its interaction with components which already exist in the MOEA/D environment.

Adding new components to MOEADr
In addition to a wide variety of components from the literature, the MOEADr package also supports the use of user-defined components. For this example, we will create a simple variation operator that does not exist in the package in its current state, and compare it with the original MOEA/D. Consider the following "Gaussian Mutation" operator: given a set X of solutions x i ∈ X we add, with probability p, a noise r ij ∼ N (µ, σ) to each x ij ∈ x i ∈ X. The R code for this operator is as follows: R> variation_gaussmut <-function(X, mean = 0, sd = 0.1, p = 0.1, ...) { R> # vector of normally distributed values, R> #length(R) = nrow(X) * ncol(X) R> R <-rnorm(length(X), mean = mean, sd = sd) R> # Apply random binary mask, probability = p R> R <-R * (runif(length(X)) <= p) R> # Add mutations to the solution matrix R> return (X + R) R> } We would like to highlight a few characteristics of the code sample above. First, the name of the function must be in the form variation [functionname]. The MOEADr package uses function name prefixes to perform some automated functions such as listing existing components and error checking. The list of current function prefixes and their meaning is: • constraint : Constraint handling components Second, the parameters in the definition of the variation operator function must include: the solution set matrix X, any specific parameters for the function, and finally an ellipsis argument to catch any other inputs passed down by the main moead() function, such as objective values and former solution sets.
Extensively commented examples of using these parameters are available in the source code for the variation operators included in the package, such as the Binomial Recombination operator (variation binrec()). Other component classes follow similar rules, as documented in the Writing Extensions for the MOEADr Package vignette.
If a given function is not available in the MOEADr package environment, it will search for it in the base R environment. Therefore, if you have named your component correctly, all you need to do is add it to the appropriate parameter in the moead() call.

Summary and conclusions
In this article we presented the R package MOEADr, implementing a new, componentbased formulation of the MOEA/D framework. This formulation breaks down the algorithm into independent components that can be separately replaced or configured. We described many recent proposals in the MOEA/D literature in an unified mathematical formulation fitting this framework.
The package allows users to easily instantiate a large variety of works from the MOEA/D literature, which facilitates the reproducibility of published re-sults and the application of known variants to problems of interest. By implementing the components listed in Table 4 and allowing the definition of new functions for any component of the algorithm, the package also allows practitioners to quickly test new proposals and ideas.
Recently, [4] introduced a component-wise framework for dominance-and indicator-based MOEAs (a vision that is partially implemented in the ecr package by [8]), and described how that framework was useful for the automated generation of algorithms. In the present work we introduced a counterpart for decomposition-based approaches, covering the third major branch of multiobjective evolutionary algorithms.
Besides expanding the library of available components, current work for improving the package includes the addition of parameter self-adaptation and of performance-based stop criteria. The possibility of including preference information within the MOEADr framework to bias the search towards specific subsets of the Pareto-optimal front [25,26], as well as expanding the package to deal with robust multiobjective optimization problems [24] are also among the improvements planned for a future release.