Decomposition-based evolutionary dynamic multiobjective optimization using a difference model

Applied Soft Computing Journal 76 (2019) 473–490

Contents lists available at ScienceDirect

Applied Soft Computing Journal journal homepage: www.elsevier.com/locate/asoc

Decomposition-based evolutionary dynamic multiobjective optimization using a difference model ∗

Leilei Cao a , Lihong Xu a , , Erik D. Goodman b , Hui Li c a

Department of Control Science and Engineering, Tongji University, Shanghai, 201804, China BEACON Center for the Study of Evolution in Action, Michigan State University, East Lansing, MI, 48824, USA c School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, 710049, China b

highlights • • • • •

Each Pareto-optimal solution was assumed to have the same motion with the centroid. Two or three historical locations of the centroid were collected. A difference model was proposed to predict the motion of the centroid. Part of old solutions were retained to mix with the predicted solutions. MOEA/D framework was modified to adapt the dynamic change.

article

info

Article history: Received 5 March 2018 Received in revised form 7 November 2018 Accepted 27 December 2018 Available online 2 January 2019 Keywords: Decomposition Difference model Evolutionary algorithm Dynamic multiobjective optimization

a b s t r a c t This paper presents a novel prediction model combined with a multiobjective evolutionary algorithm based on decomposition to solve dynamic multiobjective optimization problems. In our model, the motion of approximated Pareto-optimal solutions (POS) over time is represented by the motion of the centroid, and the other solutions are assumed to have the same motion as the centroid. A history of recent centroid locations is used to build a difference model to estimate the later motion of the centroid when an environmental change is detected, and then the new locations of the other solutions are predicted based on their current locations and the estimated motion. The predicted solutions, combined with some retained solutions, form a new population to explore the new environment, and are expected to track the new POS and/or Pareto-optimal front relatively well. The proposed algorithm is compared with four state-of-the-art dynamic multiobjective evolutionary algorithms through 20 benchmark problems with differing dynamic characteristics. The experimental studies show that the proposed algorithm is effective in dealing with dynamic problems and clearly outperforms the competitors. © 2018 Elsevier B.V. All rights reserved.

1. Introduction

dynamic changes in number of objectives [9,10]. In this paper, we only focus on the following type of DMOPs [11]:

Evolutionary dynamic multiobjective optimization (EDMO) has recently attracted intense interest in the evolutionary computation research community [1,2]. Due to the uncertainties and dynamics of many real-world multiobjective optimization problems [3], an EDMO algorithm must not only find the Pareto-optimal front (POF) and the Pareto-optimal solutions (POS) at a given time, but also track the moving POF and/or POS [4–6]. There are several types of dynamic multiobjective optimization problems (DMOPs), distinguished according to the nature of the dynamics exhibited [7,8]— e.g., those with dynamic objective functions, or dynamic constraints, and/or dynamically changing decision variables, or even

min F (x, t ) = (f1 (x, t ) , f2 (x, t ) , . . . , fm (x, t ))T

∗ Corresponding author. E-mail addresses: [email protected] (L. Cao), [email protected] (L. Xu), [email protected] (E.D. Goodman), [email protected] (H. Li). https://doi.org/10.1016/j.asoc.2018.12.031 1568-4946/© 2018 Elsevier B.V. All rights reserved.

subject to x ∈ Ω

(1)

where m is the number of objectives, t = 0, 1, 2. . . represents discrete time instants, x = (x1 , x2 , . . . , xn )T is the decision variable vector, in which n is the number of decision variables, and Ω represents the decision space. The objective vector F (x, t ) consists of m time-varying objective functions that change intermittently. Research on EDMO involves various aspects, including construction of benchmark functions and performance metrics, development of algorithms, application of algorithms to real-world problems, etc. Over the past few years, there have been a variety of contributions to these aspects [12]. Developing algorithms for dealing with DMOPs is the main focus of this research. Many algorithms on EDMO are derived from multiobjective evolutionary

474

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

algorithms (MOEAs) that were originally developed for static multiobjective optimization problems—e.g. NSGA-II [13], MOEA/D [14, 15], RE-MEDA [16], and COEA [17]. In order to adapt them to dynamic environments, these MOEAs have been modified in two respects: (1) enhancing their convergence speed and (2) introducing more diversity of individuals after detecting an environmental change [18,19]. One formulation of the task of a dynamic MOEA is to approximate the POF and POS prior to each environmental change [20], thus tracking the changing POF and/or POS as closely as possible. This requires that the algorithm converges rapidly in a static environment, which is also crucial to its tracking capability in a dynamic environment. Most dynamic MOEAs in recent literature increase the diversity of individuals via the following strategies: randomly reinitializing [21], parameter tuning [22], memory-based reinitializing [23–25], or prediction-based reinitializing [26–28]. Generally, the prediction-based strategy is associated with the memorybased strategy, because the prediction model is based on historical information about the population [29]. In many real-world DMOP formulations, the objective functions of problems change according to some regular rules, rather than randomly between two consecutive environments [12]. Also, in many cases, the problems change smoothly, rather than rapidly, in most situations. Hence, the prediction-based strategy has been more and more adopted to introduce the predicted solutions into the re-initialized population, in order to approximate the new POF and/or POS. Although there are a certain number of dynamic MOEAs using the prediction-based strategy, more studies of this strategy are greatly needed—e.g., what kinds of prediction models are proper to predict new locations of individuals, how many historical solutions are needed to build the prediction model, and how many individuals in the current population should accept the predicted new locations. In addition, the proposed prediction models in reported dynamic MOEAs appear to be few—a linear prediction model and an autoregression model are the most commonly used. Simpler yet effective prediction models should be presented to deal with dynamic changes in DMOPs. In this paper, a novel prediction model is proposed, combined with a MOEA based on decomposition (MOEA/D) to solve DMOPs. MOEA/D decomposes a multiobjective optimization problem into a number of single-objective optimization subproblems through an aggregation function, then optimizes them simultaneously under a specific environment [14]. It has received wide attention since its introduction. In this paper, several sets of solutions obtained by MOEA/D under previous environments are used to build the prediction model. Motivated by the idea of the centroid [27,30], a sequence of solutions’ centroids is used to estimate the next motion through a difference model that is invoked whenever an environment change is detected, to estimate the moving trend of solutions. The other members are assumed to have similar moving trends with the centroid. Therefore, their new locations can be predicted at least approximately by their current locations relative to the centroid and the estimated motion of the centroid. The predicted solutions are mixed with some solutions that were retained from the most recent environment to form a new population, and exploration starts in the new environment. The major contributions of this paper include the following: (1) A novel prediction strategy that is based on a difference model has been presented, in which a first-order difference model has been extended to a second-order difference model to predict the moving trends of centroid of the solutions obtained. (2) The implementation of the novel prediction strategy combined with MOEA/D has been suggested, including which individual(s) and how many individuals should be predicted. (3) Experiments have been conducted to compare the proposed algorithm with four state-of-the-art dynamic MOEAs, including analyzing their performances on several benchmark functions.

The remainder of this paper is organized as follows. Section 2 presents background of this work, including a brief overview of some EDMO algorithms in recent literature. Section 3 presents our proposed algorithm, including the novel prediction model and its combination with MOEA/D. Experimental design based on a set of benchmark functions is given in Section 4. A comparative study and further analysis are presented in Section 5. Finally, the discussion and conclusion are given in Section 6. 2. Background 2.1. Challenges of EDMO Changes in the objective functions of an MOP introduce two novel situations for MOEAs. The first is that Pareto-optimal solutions approximated in a new environment may be dominated by previous, no longer optimal, solutions remaining from an earlier environment, even though the new solutions are, in fact, Pareto optimal in the new environment, if the Pareto-optimal front also changes. This could lead to improper discarding of members of the POS in the new environment, until the environmental change is recognized and the status of previously POS members is changed. The second situation is that a POS that has been stable may suddenly have some or all of its solutions dominated by newly found solutions, if the POF has also changed. This situation is potentially a useful indicator of an environmental change, meaning that the non-dominated status of any previous solutions still retained should now be re-evaluated. The domination of POS members, given that they have already converged adequately, can serve to signal the MOEA that an environmental change has occurred, and that it should seek a new POS as soon as possible, disregarding the domination status of previously ‘‘known’’ POS members, which may no longer be Pareto optimal. Unlike what is frequently seen in solving dynamic singleobjective optimization problems, the individuals in the population may not completely lose diversity, which means that the current individuals can still generate offspring and explore the search space, so long as they are no longer assumed to be Pareto-optimal. To illustrate the challenges of EDMO clearly, we first use FDA1 [31] as the benchmark problem and MOEA/D-DE [32] as the multiobjective optimization algorithm for investigation.1 In FDA1, the POS changes over time while the POF does not change. For initial simplicity, we assume that the moment of changes are known, and all individuals in the population are re-evaluated and the reference point is updated at the beginning of the next environment. We only investigate whether the current individuals could continue to explore the search space; therefore, no randomly reinitialized solutions are introduced into the population. Fig. 1(a) shows the population distribution in the objective space obtained by MOEA/D at the end of the first time window (just before the change), and Fig. 1(b) shows the population distribution obtained at the end of the second time window. Obviously, the algorithm cannot well approximate the true POF within the given number of evaluations before the change. The population distribution in the objective space illustrates that the individuals still maintain a certain degree of diversity, and that they can subsequently generate diverse offspring and continue to find the new POS. At the end of the second time window, the algorithm obtains a better approximation to the true POF than that found in the first time window. This illustrates that the individuals can search for the new POS in the second time window based on the former solutions 1 Here the severity of change and the frequency of change in FDA1 are both set to 10, and the dimension of decision variables in FDA1 is set to 10. In MOEA/D-DE, CR = 0.5, and F = 0.5; The population size is set to 100, and for the other parameters, refer to [32].

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

475

Fig. 1. Population distribution in the objective space obtained by MOEA/D on FDA1 at the end of time window.

found just before the change. The former solutions are helpful in searching for the new POS in the later time window. However, the algorithm still cannot well approximate the true POF within the given number of evaluations, even though the initial population in the second time window derives from the former solutions in the previous time window. The above example illustrates that, for at least some EDMO problems, increasing diversity of individuals is not critical to tracking the changing POS; rather, leading the individuals to approximate the changing POS and/or POF as quickly as possible is more important. That is the biggest challenge posed by EDMO. In addition, for problems with intermittent changes, detecting the moment of change is another challenge for the algorithm. We focus only on the first challenge in this paper. 2.2. Prediction for EDMO Many approaches have been proposed in order to lead the individuals to track the moving POF and/or POS, e.g., partially randomly reinitializing the population, tuning the parameters, memory-based strategies [33,34], and prediction-based strategies. In this paper, we mainly investigate the prediction-based strategies for EDMO. A prediction-based strategy is based on an assumption that the problems change at least somewhat regularly, according to some rules, rather than changing with no pattern between successive environments. Otherwise, a prediction-based strategy may not be more effective than other strategies when dealing with DMOPs. Hatzakis and Wallace [28] proposed a feed-forward prediction strategy, in which an autoregressive model was used as a forecasting method. Zhou et al. [27] also chose the univariate autoregression model as the prediction method. In their algorithm, the POS set was divided into two parts: a center point and a manifold. The prediction model was built through a series of historical centers to estimate the new center in the next time window. Koo et al. [30] used a weighted average approach based on the history of previously discovered solutions to estimate the predicted direction and magnitude of the next change. Liu et al. [35] used a linear regression prediction strategy to predict most of the individuals in the population at the new moment, and the remaining individuals were randomly reinitialized to increase the diversity. Wu et al. [36] proposed a directed search strategy to improve the performance of NSGA-II in dynamic environments. In their algorithm, the idea of center points was adopted to estimate the direction of motion of Pareto-optimal solutions by a linear vector difference, which was used to generate new individuals that were also perturbed by Gaussian noise. This idea was also applied by Jiang and Yang [12] in their proposed algorithm. Muruganantham et al. [11,37] employed

a Kalman filter (KF) combined with MOEA/D to solve DMOPs. They designed two variants of KFs for prediction, a 2-D KF and a 3-D KF. The 2-D KF model was a first-order linear model perturbed by Gaussian noise, and the 3-D KF model was a second-order linear model perturbed by Gaussian noise. Both variants are discrete and linear KFs, which were used to predict the locations of the POS after a change is detected. 2.3. Weakness of existing prediction-based strategies Most of these reported prediction-based strategies have something in common, in that the prediction model is perturbed by Gaussian noise. The reason for involving Gaussian noise is the existence of prediction error. Generally, there are two sources of error in the prediction model [29]: input and output. (1) The input of the prediction model is a time series of the approximated POS in the previous environments. However, these approximated POS may deviate from the path under the unknown governing dynamics because the algorithm has not converged properly (the dynamics have not been well approximated). Therefore, the deviation causes input error. (2) Even if the input data is accurate, the prediction model may output an inaccurate result. The prediction model is built through a time series of the approximated POS, for which it is hard to estimate all types of changes. For example, in [36], the essential principle of their prediction model is a linear model, which is based on an assumption that the problems change over time in a linear pattern. However, the problems may change in other patterns; thus the predicted location may deviate from the true location in the next environment. Although prediction error certainly exists, these reported approaches may not take full advantage of the effectiveness of MOEAs. The main aim of a prediction model in this context is to provide an initial population for the MOEAs, and they need not be in the exact locations of the new POS. As long as the predicted locations are close to the true locations of the new POS in the next environment, and then the MOEAs can explore the search space on the basis of the predicted locations. However, the involvement of Gaussian noise may cause the predicted location to be far away from the true location, which causes the MOEAs to spend more evaluations on discovering the new locations of the POS in the new environment. In addition, the input error may not exhibit a Gaussian distribution and it is difficult to build an accurate model to describe the input error. Consequently, unlike the other reported prediction-based strategies, the involvement of Gaussian noise will be abandoned in our prediction model. Besides involving Gaussian noise, most of these reported prediction-based strategies ignore retaining the approximated POS

476

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

just before the change when updating the population after the change. Sometimes, the change severity of DMOPs is relatively slight, which means the new locations of the true POS may be close to the previous locations just before the change. In that way, the approximated POS in the former environment may be useful to the new population for tracking the new locations of the POS. Accordingly, the former approximated POS just before the change will be retained to update the population in our algorithm. 3. Proposed algorithm In order to predict the locations of the POS, we have to make some assumptions, as are made in [27]. For the DMOP’s to be considered here: (1) the problem remains stationary during a fixed time span; (2) the POS set or POF of consecutive spans change according to some simple rules. Based on these two assumptions, we propose a decomposition-based evolutionary dynamic multiobjective optimization using a difference model (MOEA/D-DM) to solve such DMOPs. The main steps of the proposed MOEA/D-DM are described as follows. Step 0. Initialize a population P, and initialize the reference point z∗. Step 1. Check for an environmental change. Step 1.1. If the environment changes, go to Step 1.2, else go to Step 2. Step 1.2. Predict the next location of selected individuals, and update the population. Step 1.3. Re-evaluate the new population, and update the reference point z ∗ . Step 2. Perform reproduction and updating of the population in this environment. Step 3. If the stopping criterion is satisfied, stop; else go to Step 1. The main steps of MOEA/D-DM are based on the framework of MOEA/D-DE [32], which is a popular and typical multiobjective evolutionary algorithm, in which a differential evolution (DE) operator and a polynomial mutation operator are used for producing new solutions. To the base of MOEA/D-DE, Step 1 is added— i.e., checking for an environmental change and, if one is found, predicting the new locations of POS, and updating the new population. This paper mainly focuses on Steps 1.2 and 1.3, and for the details of the other steps, the reader is referred to [32]. In many real-world DMOPs, we have no prior knowledge of the environmental change, including the moment of change, or the severity of change. In order to detect the moment of change, some methods have been proposed elsewhere—e.g., re-evaluating a fraction of the existing solutions during each generation [21,38,39], or assessing some statistical information about some selected population members [31]. Although the first method is easy to implement, it is based on an assumption that there is no noise in the function evaluations. The second method needs some additional problemdependent parameters. In this paper, we use the first method to detect the environmental change, which is implemented at the beginning of each generation. If their average fitness value in each objective changes, we conclude that the environment has changed. In the following subsections, we focus on how to predict the new locations of the POS and how to update the new population after detecting the environmental change. 3.1. MOEA/D MOEA/D decomposes the problem (1) at any time instant into a number of single-objective optimization subproblems through an aggregation function [40,41]. In this paper, the Tchebycheff approach acts as the decomposition method. Let λ1 , . . . , λN be a set of even spread weight vectors, a multiobjective optimization

problem can be decomposed into N scalar optimization problems, and the jth subproblem is j

minimize g x ⏐λj , z ∗ = max λi ⏐fi (x) − zi∗ ⏐

( ⏐

)

⏐

⏐

1≤i≤m

subject to x ∈ Ω

(2)

where z = (z1 , . . . , zm ) is the reference point, i.e., zi = min{fi (x) , x ∈ Ω } (for a minimization problem) for each i = 1, . . . , m. MOEA/D minimizes these N subproblems simultaneously in a single run. j In MOEA/D, the neighborhood of weight vector { 1λ is defined } as the set of its several closest weight vectors in λ , . . . , λN . The neighborhood of the jth subproblem consists of all the subproblems with the weight vectors from the neighborhood of λj . ∗

∗

∗ T

∗

3.2. Difference model-based prediction The difference model is built through a series of POS obtained in several previous time windows (also called environments). We assume that a series of POS obtained in the previous time windows are denoted as PT −K +1 , . . . , PT , and that the POS for time window T + 1 can be considered as a function of PT −K +1 , . . . , PT [26]: PT +1 = fun(PT −K +1 , . . . , PT , t)

(3)

where PT +1 denotes the POS for time window T + 1, and K represents the number of previous time windows that PT +1 depends on in the prediction model. We assume that each set of POS (PT −K +1 , . . . , PT ) is composed of N solutions; then xti (i = 1, 2 . . . , N , t = T − K + 1, . . . , T ) represents the ith solution in Pt , t = T − K + 1, . . . , T . In MOEA/D, each subproblem that is determined by a particular weight vector has a corresponding solution, which means that xti (i = 1, 2 . . . , N , t = T − K + 1, . . . , T ) represents the solution of the ith subproblem of a MOP in the tth environment. Accordingly, xiT −K +1 , . . . , xTi are a series of solutions that describes the motion of the solution in the ith subproblem over time. A prediction model can estimate its next solution in the following time window from the former solutions. Theoretically, each solution should be estimated by a separate prediction model; therefore, there will be N models to describe the motion of N solutions in the POS exactly. However, it is too complex to use so many different prediction models. It is assumed that the N solutions have similar dynamic characteristics, and the aim of the prediction model is only to provide the initial population for the evolutionary algorithm; thus a generic prediction model is sufficient to describe the motion of these solutions. Although each new solution at time window T + 1 can be predicted through its corresponding former solutions, how can we guarantee the accuracy of these historical solutions? As mentioned in Section 2.3, there exists input error. In addition, each solution that is obtained during a series of evaluations may deviate in its trajectory over time windows in different ways. Thus, inaccurate historical solutions may produce a predicted result far from the actual one. To reduce the interference from the input error, the centroid of the population is used to describe the motion of the solutions over time [12,27,30,36]. We assume that each member in the population has the same motion trend over time as that of the centroid. However, this condition is too strict; all solutions may have similar motions to the centroid, but not exactly the same in most situations in nature. However, given that the predicted new positions are only starting points for the next search process, the prediction error can often be tolerated, as demonstrated below. Let C T be the centroid in PT , then C T can be computed as: CT =

1 ∑

|PT |

x

x∈PT

where |PT | is the cardinality of PT , and x is a solution in PT .

(4)

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

477

real-world optimization problems may change in various other simple ways—e.g., nonlinear or stochastic [28]. Thus the secondorder difference model considers a varying motion for the centroid in different time windows. Fig. 3(b) describes the predicted positions of the centroid by the second-order difference model in a 2-

−−−−−−−−→ → → D decision space. − v1 represents C T −1 − C T −2 , and − v2 represents

−− −−−−→ C T − C T −1 . The second-order difference model considers the centroid’s motion as a ‘‘uniformly accelerated motion’’ over discrete time. After predicting the moving trend of centroid C T at the end of time window T, the potential locations of other members in PT at time window T + 1 can be predicted as follows:

Fig. 2. Motion of solutions at different time windows in a 2-D decision space.

−−→

xTi +1 = xTi + ∆X T , i = 1, 2 . . . , N Fig. 2 illustrates the posited motions of solutions at various time windows in a 2-D decision space. As mentioned above, each solution is assumed to have similar motion to that of the centroid. Thus the motions of solutions are described by those of the centroids C T −3 , C T −2 , C T −1 , C T in Fig. 2. The prediction model is

−−→

designed to estimate the moving trend ∆X T of C T at the end of time window T through positions of the centroid at several previous time windows, in order to estimate the potential positions of other solutions at time window T + 1. The range of previous time

−−→

windows that the moving trend ∆X T depends on can start from 1 to T, theoretically. However, involving too many historical centroids may cost excessive computational resources and could be misleading, if the process is changing non-uniformly. Considering the tradeoff between computational resource requirements and prediction accuracy, we have designed two schemes to predict the

−−→ −−→T that ∆X depends on:

moving trend ∆X T , according to the number of previous centroids

−−→

−−−−−−→

First-order difference: ∆X1T = C T − C T −1

−−→

−−−−−−→

(5)

−−−−−−→

Second-order difference: ∆X2T = C T − C T −1 + (C T − C T −1

−−−−−−−−→ − C T −1 − C T −2 )

(6)

The first-order difference model is based on an assumption that the centroid’s moving trend at the end of time window T is the same as the centroid’s motion from time window T − 1 to T, including the same step size and the same direction (denoted by a vector). Under this model, the centroid’s motion can be considered as ‘‘uniform motion’’ in different time windows, based on the assumption that the optimization problem changes linearly over time. Fig. 3(a) describes the predicted positions of the centroid by the first-order difference model in a 2-D decision space. However,

(7)

After predicting new solutions, the population should be updated in order to adapt to the new environment. 3.3. Update of the new population Unlike most EDMO algorithms in the literature, some solutions obtained in the latest former environment are retained in the new population in our algorithm. In real-world DMOPs, the POS of consecutive MOPs are often similar to each other [7]. Therefore, these old solutions found just before the change can also be added to the initial population in the new environment. The old solutions may perform better than an entirely reinitialized population as is used in most current DMOEAs. The new population is composed of two kinds of solutions: the old solutions and the predicted solutions. In MOEA/D, all solutions found just before the change have been automatically arranged according to the uniformly distributed weight vectors with which they are associated. Since each subproblem is optimized by using information from its neighboring subproblems [42], the old solutions and the predicted solutions both tend to be distributed somewhat uniformly. We use a parameter q (q ∈ [0, 1]) to control the number of predicted solutions in the new population. Thus there will be q ∗ N solutions that are generated by accepting the prediction results, and the others are retained from the solutions obtained just before the change. Fig. 4 shows an example of the distribution of a new population in a two-objective space at the beginning of a new time window. In Fig. 4, the neighborhood of each subproblem consists of both the old solution and the predicted solution. The resulting collaboration between the old solution and the predicted solution tends to lead the population to track the new POS and/or POF relatively rapidly. If a particular neighborhood includes only old solutions, then the prediction loses its effect for that neighborhood, at least in the near term.

Fig. 3. Positions of centroid predicted by different models in a 2-D decision space.

478

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Fig. 4. Distributions of a new population in a two-objective space at the beginning of a new time window.

tuned to possess a number of challenging characteristics, including a mixed Pareto-optimal front, nonmonotonic and time-varying variable linkages, and mixed types of changes. The time instance t involved in these problems is defined as t = 1/nt ∗ ⌊τ /τt ⌋ (where nt , τt , and τ represent the severity of change, the frequency of change, and the iteration counter, respectively). The definitions of these problems can be found in Tables 1–3. There are a number of metrics to be used for performance assessment of DMOEAs—e.g., Inverted Generational Distance (IGD) [44,45] , Schott’s Spacing Metric [46], Maximum Spread [47] and Hypervolume Difference [48]. In this paper, a modified version of the IGD (MIGD) metric is used to assess the algorithms, as is suggested in [11] and [27], given that the true Pareto front sequences of these benchmarks are known. ∗ Let Qt be a set of uniformly distributed points in the true POFt , and Qt be an approximation of POFt . The IGD metric is defined as:

∑ IGD Q t ∗ , Q t =

(

)

v∈Q t ∗

d(v, Q t )

|Q t ∗ |

(8)

where d v, Q t ⏐ = min ⏐ u∈Q t ∥F (v) − F (u)∥ is∗ the distance between v and Qt , and ⏐Q t ∗ ⏐ is the cardinality of Qt . The IGD metric can measure both diversity and convergence. To have a low IGD value, Qt should be very close to the true POFt and cannot miss any part of the whole POFt . In the experiments, 500 points and 1000 points are uniformly sampled from the true POFt of biobjective problems and triobjective problems, respectively, for computing the IGD metrics. The MIGD metric is defined as the average of the IGD values across multiple time windows over a single run:

(

MIGD =

The procedure for updating the population is given in Algorithm 1. In this procedure, 50% of the population members (an initial value for q) accept predicted solutions, and they are assigned as every other member (i.e., old and predicted solutions alternate). The prediction scheme works beginning from the third time window, in which the first-order difference model is applied to estimate the new solutions during the third time window. Starting from the fourth time window, there are three former centroids of the POS to be used. Therefore, the second-order difference model is adopted from that point forward, to better predict the new locations of the POS in the later time windows. Note that the prediction model is not available in the second time window, due to limited historical information. Consequently, all solutions are retained to track the new POS in the second time window, we avoid using a random re-initialization method for the POS. The whole pseudocode of MOEA/D-DM is given in Algorithm 2. 4. Experimental design

)

1 ∑

|T |

IGD(Q t ∗ , Q t )

(9)

t ∈T

where T is a set of discrete time instances in a run and |T | is the cardinality of T. 4.2. Compared algorithms and parameter settings Four popular DMOEAs are used for comparison in our empirical studies, including MOEA/D [32], MOEA/D-KF [11], DNSGA-II-A [21] and DDS [36]. The original MOEA/D is modified to deal with the dynamic problems. A fraction of its population is randomly reinitialized when the environment changes, which is similar to the idea of DNSGA-II-A. MOEA/D-KF used a Kalman filter prediction model to guide the search toward the changed optima, and was also based on the MOEA/D algorithm. DDS incorporated a linear prediction model along with directed local search into NSGA-II. These four algorithms that were compared also use the DE operator and the polynomial mutation operator to generate new individuals. The problems and algorithm parameter settings are as follows.

4.1. Benchmark problems and performance metric The proposed algorithm is tested on 20 benchmark problems, including five FDA problems [31], three dMOP problems [17], four ZJZ problems (F5-F8) [27], and eight JY problems (JY1-JY8) [43]. The FDA test suite is commonly used to evaluate the performance of DMOEAs. Note that FDA2 has been modified, as it was in [21]. The dMOP problems are an extension of the FDA test suite. In the FDA and dMOP test suites, the POF and/or POS change over time, and there is a linear correlation between the decision variables. In the ZJZ problems, the correlation between the decision variables is nonlinear, which challenges DMOEAs. The JY test suite is a recently proposed benchmark framework, which is able to be

4.2.1. The problem parameters To study the impact of change frequency and change severity, a variety of parameters are used. The severity of change is nt = 10, 5. The frequency of change is τt = 5, 10, 20. Therefore, each problem has 6 different cases of changes. 4.2.2. Common parameters in all algorithms The population size is set as 100 and 300 for biobjective problems and triobjective problems, respectively. CR = 0.5, F = 0.5 in the DE operator. η = 20, pm = 1/n in the polynomial mutation operator (n is the number of variables).

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

4.2.3. Parameters in MOEA/D To evaluate the MOEA/D-based algorithms fairly, some parameters in MOEA/D-DM, MOEA/D and MOEA/D-KF are set to the same values: neighborhood size is set as 20, δ = 0.8, nr = |E |. E is the mating pool, from which individuals are selected to produce a new solution. For details about E, please see [32]. Note that the value of nr in this paper is unusual, and is a key parameter in solving DMOPs; we will discuss it below. 4.2.4. Other parameters In MOEA/D and DNSGA-II-A, 20% of population members are randomly reinitialized within the domain when the environment changes. For all algorithms, a maximum of 10% of population members are chosen for re-evaluation to detect environmental

479

changes. Other parameters of DSS are as shown in [36]. The number of total generations in each run is fixed to be 40 ∗ τt , which ensures that 40 environmental changes happen in each experiment. Each algorithm is run independently 30 times for each test instance. 5. Experimental results 5.1. Results on FDA problems Table 4 shows the mean MIGD values obtained and their standard deviations over 30 runs by five algorithms on the FDA test instances. The best values are bolded. Wilcoxon’s rank sum test is used to indicate significance between results of MOEA/D-DM and

480

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Table 1 FDA and dMOP test suite. Instance

Domain

n

Remarks

[0, 1] × [−1, 1]n−1

10

POS changes POF is static

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

FDA3

f1 (x) = (x1 + x2√ )/2 f2 (x) = g(x)(1 − f∑ 1 /g(x)) n g (x) = 1 + G (t ) + i=3 (xi − G(t))2 G (t ) = |sin(0.5π t)| , F (t) = 102 sin(0.5π t) POS(t) : 0 ≤ x1 , x2 ≤ 1, xi = G (t ) , i = 3, . . . , n

[0, 1]n

10

POS changes POF changes

FDA4

f1 (x) = ((1 + g (x)) cos (0.5π x1 ) cos(0.5π x2 )) f2 (x) = ((1 + g (x)) cos (0.5π x1 ) sin(0.5π x2 )) f3 (x) = ∑ ((1 + g (x)) sin (0.5π x1 )) n 2 g (x) = i=3 (xi − G(t)) G (t ) = |sin(0.5π t)| POS(t) : 0 ≤ x1 , x2 ≤ 1, xi = G (t ) , i = 3, . . . , n

[0, 1]n

10

POS changes POF is static

[0, 1]n

10

POS changes POF changes

dMOP1

f1 (x) = x1 f2 (x) = g(x)(1 − (f1 /g(x))H(t) ) ∑n 2 9 g (x) = 1 + n− i=2 xi 1 H (t ) = 0.75∗ sin (0.5π t) + 1.25 POS(t) : 0 ≤ x1 ≤ 1, xi = 0, i = 2, . . . , n

[0, 1]n

10

POS is static POF changes

dMOP2

f1 (x) = x1 H(t) f2 (x) = g(x)(1 ∑n− (f1 /g(x)) ) g (x) = 1 + i=2 (xi − G(t))2 H (t ) = 0.75∗ sin (0.5π t) + 1.25, G (t) = |sin(0.5π t)| POS(t) : 0 ≤ x1 ≤ 1, xi = G(t), i = 2, . . . , n

[0, 1]n

10

POS changes POF changes

dMOP3

f1 (x) = xr √ f2 (x) = g(x)(1 ∑ − f1 /g) g (x) = 1 + x ∈x\xr (xi − G(t))2 i G (t ) = |sin(0.5π t)| , r = U(1, . . . , n) POS(t) : 0 ≤ xr ≤ 1, xi ∈ x\xr = G(t), i = 1, . . . , n

[0, 1]n

10

POS changes POF is static

FDA1

Description f1 (x) = x1 √ f2 (x) = g(x)(1 ∑n− f1 /g(x)) 2 g (x) = 1 + i=2 (xi − G(t)) G (t) = sin(0.5π t) POS(t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n f1 (x) = x1 (H (t )+

FDA2

∑n (x −H(j)/4)2 ) j=n−4 j

) f2 (x) = g(x)(1 − (f1 /g(x))2 ∑n−5 g (x) = 1 + i=2 x2i H (t ) = 2sin(0.5π (t − 1)) POS(t) : 0 ≤ x1 ≤ 1, xi = 0, i = 2, . . . , n − 5 xj = H (t ) , j = n − 4, . . . , n F (t )

F (t)

(

F (t)

)

cos(0.5π x2 ))

(

F (t) x1

)

F (t) x2 ))

(

F (t) x1 2

f1 (x) = ((1 + g (x)) cos 0.5π x1 f2 (x) = ((1 + g (x)) cos 0.5π FDA5

f3 (x) = ((1 + g (x)) sin 0.5π

∑n

)

F (t)

sin(0.5π

)

g (x) = G (t ) + i=3 (xi − G(t)) G (t ) = |sin(0.5π t)| , F (t) = 1 + 100 sin4 (0.5π t) POS(t) : 0 ≤ x1 , x2 ≤ 1, xi = G (t ) , i = 3, . . . , n

each competitor at the 5% significance level. The results of statistical tests are marked as ‘‘+’’, ‘‘−’’, or ‘‘∼’’ according to whether MOEA/D-DM is statistically better than, worse than, or not statistically different from the corresponding algorithm, respectively. It is obvious that MOEA/D-DM performed best on the majority of the FDA instances, implying that it has the best ability to track changing POS and/or POF in most cases. The results show that the three MOEA/D-based algorithms performed clearly better than DNSGA-II-A and DDS on FDA1-FDA4, which illustrates that the MOEA/D-based algorithms have better tracking abilities than the NSGA-II-based algorithms on these problems. The performances of all algorithms were closely related to the change frequency and change severity. The fast and high-severity change (τt = 5, nt = 5) challenged all algorithms, due to the limited number of evaluations allowed and the large shift between two consecutive environments for a solution in the decision space. When more evaluations are allowed (τt = 10 and 20), all algorithms obtain better results. In addition, they perform better when the change severity is relative slight. FDA1 is a simple problem, in which the POS changes over time as a sine pattern but the POF remains stationary. In FDA2, the POF changes from convex to non-convex shapes in the objective

space over time and a subset of the decision variables also changes. MOEA/D-DM has obtained quite good results when τt = 10 and 20 for these two problems. MOEA/D-KF performs better than MOEA/D on FDA1 and DDS also performs better than DNSGA-II-A on this problem, which indicates that a prediction scheme can effectively enhance the search efficiency of the corresponding algorithm on this problem. The approach of partial random re-initialization is successful for MOEA/D and DNSGA-II-A on FDA2, which illustrates that the partial random re-initialization is also even more effective than some prediction schemes on some problems. FDA3 is a problem in which the environmental changes shift the POS and affect the density of points on the POF. Performances of all algorithms on this problem deteriorate comparing with their respective performances on FDA1 and FDA2. FDA4 and FDA5 are dynamic triobjective optimization problems; in FDA4 only POS changes over time, whereas both POS and POF change over time in FDA5. MOEA/DDM still has slight superiority compared to MOEA/D and MOEA/DKF on FDA4. Note that DDS performs worse than DNSGA-II-A on this problem. All algorithms seem less sensitive to the severity of change on FDA4. No single one of the algorithms beat all others on all cases of FDA5. Since the density of the POF changes over

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

481

Table 2 ZJZ test suite. Instance

Description f1 (x) = |x1 − a| + H

∑

2 i∈I1 yi

f2 (x) = |x1 − a − 1|H +

∑

i∈I2

Domain

n

Remarks

[0, 5]n

10

POS changes POF changes

[0, 5]n

10

POS changes POF changes

[0, 5]n

10

POS changes POF changes

[0, 1]2 × [−1, 2]n−2

10

POS changes POF changes

y2i

i n

F5

yi = xi − b − 1 + |x1 − a|H + H = 0.75 sin (π t ) + 1.25 a = 2 cos (π t) + 2, b = 2 sin (π t) + 2 I1 = {i |1 ≤ i ≤ n, i is odd } , I2 = {i |1 ≤ i ≤ n,(i is ev en } ) POS(t) : a ≤ x1 ≤ a + 1, xi = b + 1 − |x1 − a|

H + ni

, i = 2, . . . , n

f1 (x) = |x1 − a|H + i∈I y2i 1 ∑ f2 (x) = |x1 − a − 1|H + i∈I y2i 2

∑

i

F6

yi = xi − b − 1 + |x1 − a|H + n H = 0.75 sin (π t ) + 1.25 a = 2 cos (1.5π t) sin (0.5π t) + 2 b = 2 cos (1.5π t) sin (0.5π t) + 2 I1 = {i |1 ≤ i ≤ n, i is odd } , I2 = {i |1 ≤ i ≤ n,(i is ev en } ) POS(t) : a ≤ x1 ≤ a + 1, xi = b + 1 − |x1 − a|

f1 (x) = |x1 − a| + H

∑

2 i∈I1 yi

f2 (x) = |x1 − a − 1|H +

∑

i∈I2

H + ni

, i = 2, . . . , n

y2i

i n

F7

yi = xi − b − 1 + |x1 − a|H + H = 0.75 sin (π t ) + 1.25 a = 1.7(1 − sin (π t)) sin (π t) + 3.4 b = 1.4(1 − sin (π t)) cos(π t) + 2.1 I1 = {i |1 ≤ i ≤ n, i is odd } , I2 = {i |1 ≤ i ≤ n,(i is ev en } ) POS(t) : a ≤ x1 ≤ a + 1, xi = b + 1 − |x1 − a|

F8

H + ni

, i = 2, . . . , n

f1 (x) = (1 + g(x)) cos(0.5π x1 ) cos(0.5π x2 ) f2 (x) = (1 + g(x)) cos(0.5π x1 ) sin(0.5π x2 ) f3 (x) = (1 + g(x)) sin(0.5π x1 )

∑n

(x

+x

)H

1 2 − G)2 g ( x) = i=3 (xi − 2 G = sin (0.5π t) , H = 0.75 sin (π t) + 1.25

POS(t) : 0 ≤ x1 , x2 ≤ 1, xi =

( x1 +x2 )H 2

+ G, i = 3, . . . , n

time, a good distribution of solutions on the POF found in the former environment will no longer be a good one after a change. No algorithm had evident superiority on this problem across different frequencies and severities of changes. Fig. 5(a)–(e) present the tracking of the IGD values with the environmental change obtained by the five algorithms on the FDA test instances with nt = 10, τt = 10, showing the average values over 30 runs. Observing from Fig. 5(a), (b) and (d), the tracking ability of MOEA/D-DM was most stable with environmental changes across different time windows, and MOEA/D-DM had the best IGD values at most time windows on these three problems. The IGD values obtained by DDS and DNSGA-II-A fluctuate very widely on FDA1. MOEA/D performs similarly to MOEA/D-KF on FDA1 in terms of their tracking stabilities. However, the IGD values obtained by MOEA/D-KF fluctuate a lot on FDA2 at some time windows, as do those of DNSGA-II-A. None of the algorithms can stably track the changed POS and POF on FDA3, as can be seen in Fig. 5(c), in which all IGD curves fluctuate widely. The IGD values achieved by DDS on FDA4 are obviously worse than those of the others across all changes. The three MOEA/D-based algorithms have similar performances on this problem. Similar to the performances of all algorithms on FDA3, all algorithms perform unstably on FDA5, looking across all changes. The tracking curves of these algorithms on FDA test instances indicate that the problems with changing POS and changing POF challenge DMOEAs, especially those problems with a varying density of solutions on the POF over time (like FDA3 and FDA5). 5.2. Results on dMOP The average MIGD values, standard deviations and statistical test results obtained by five algorithms on dMOP test instances are presented in Table 5. Only the POF is dynamic in dMOP1, whereas

both the POS and POF change over time in dMOP2. MOEA/D-DM performs best on all cases of these two problems. MOEA/D and MOEA/D-KF also perform well compared to DNSGA-II-A and DDS on dMOP1 and dMOP2. Like FDA1, only the POS changes in dMOP3, but the difference is that the variable in dMOP3 that controls the spread of the POF changes randomly over time. The randomness of change challenges the prediction scheme. The partial random re-initialization applied in MOEA/D shows its superiority when the environment changes quickly (τt = 5) on this problem. The Gaussian noise in the Kalman Filter model of MOEA/D-KF produces its effect when the environment changes relatively slow (τt = 10 and 20) on dMOP3. MOEA/D-DM’s performance is not statistically different from that of MOEA/D-KF when τt = 20. The results obtained by the three MOEA/D-based algorithms indicate that the prediction scheme does not show its superiority if the environment changes quickly and there are random changes in the decision variables. However, the prediction scheme with associated noise has better performance if the environment changes slowly. Note that DNSGA-II-A also performs better than DDS on dMOP3 with different frequencies and severities of changes. Fig. 5(f)–(h) present the tracking curves of these five algorithms on dMOP test instances with nt = 10, τt = 10. MOEA/D-DM and MOEA/D have similar tracking abilities on dMOP1, as seen in Fig. 5(f). Both of them respond to changes stably and have low IGD values after the third time window, which illustrates that these two algorithms have quite good convergence performances and tracking abilities on this problem. Since dMOP2 has varying POS and POF, none of the algorithms can respond stably. However, MOEA/D-DM still has relatively low IGD values in most time windows. Although no algorithms respond stably on dMOP3, the tracking abilities of MOEA/D and MOEA/D-KF are better than those of the other algorithms.

482

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Table 3 JY test suite. Instance

Description

Domain

n

Remarks

JY1

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 )) f2 (x) = ∑ (1 + g(x))(1 − x1 + Asin(W π x1 )) n 2 g (x) = i=2 (xi − G(t)) , G(t) = sin(0.5π t) A = 0.05, W = 6 POS (t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF is static

JY2

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 )) f2 (x) = ∑ (1 + g(x))(1 − x1 + Asin(W π x1 )) n 2 g (x) = i=2 (xi − G(t)) , G(t) = sin(0.5π t) A = 0.05, W = ⌊6 sin(0.5π (t − 1))⌋ POS (t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

JY3

f1 (x) = (1 + g(x))(y1 + Asin(W π y1 )) f2 (x) = ∑ (1 + g(x))(1 − y1 + Asin(W π y1 )) n 2 2 g (x) = .05 i=2 (yi − yi−1 ) , A = 0⌊ ⌋ W = ⌊6 sin(0.5π (t − 1))⌋ , α = 100 sin2 (0.5π t) y1 = |x1 sin(2α + 0.5)π x1 | , yi = xi , i = 2, . . . , n √ POS (t) : 0 ≤ x1 ≤ 1, xi = xi−1 , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

JY4

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 )) f2 (x) = ∑ (1 + g(x))(1 − x1 + Asin(W π x1 )) n 2 g (x) = i=2 (xi − G(t)) , G(t) = sin(0.5π t) A = 0.05, W = 101+|G(t)| POS (t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

JY5

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 )) f2 (x) = ∑ (1 + g(x))(1 − x1 + Asin(W π x1 )) n 2 g (x) = i=2 xi , A = 0.3 sin(0.5π (t − 1)), W = 1 POS (t) : 0 ≤ x1 ≤ 1, xi = 0, i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS is static POF changes

JY6

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 )) f2 (x) =∑ (1 + g(x))(1 − x1 + Asin(W π x1 )) n 2 g(x) = i=2 (4yi − cos (K π yi ) + 1) A = 0.1, W = 3, K = 2∗ ⌊10 ∗ |G(t)|⌋ G (t) = sin (0.5π t) , yi = xi − G (t ) , i = 2, . . . , n POS (t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

JY7

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 ))α f2 (x) =∑ (1 + g(x))(1 − x1 + Asin(W π x1 ))β n 2 g(x) = i=2 (yi − 10 cos (2π yi ) + 10) A = 0.1, W = 3, α = β = 0.2 + 2.8∗ |G(t)| G (t) = sin (0.5π t) , yi = xi − G (t ) , i = 2, . . . , n POS (t) : 0 ≤ x1 ≤ 1, xi = G (t ) , i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS changes POF changes

JY8

f1 (x) = (1 + g(x))(x1 + Asin(W π x1 ))α f2 (x) = ∑ (1 + g(x))(1 − x1 + Asin(W π x1 ))β n 2 g (x) = i=2 xi , G (t ) = sin(0.5π t) A = 0.05, W = 6, α = β2 , β = 10 − 9.8∗ |G(t)| POS (t) : 0 ≤ x1 ≤ 1, xi = 0, i = 2, . . . , n

[0, 1] × [−1, 1]n−1

10

POS is static POF changes

5.3. Results on ZJZ The ZJZ test instances (F5-F8) have a strong nonlinear correlation between decision variables; in addition, the POS and POF change over time with strong nonlinear characteristics in the decision space and objective space, respectively. Therefore, the ZJZ problems greatly challenge DMOEAs, especially when the problem changes frequently and severely. Table 6 presents the results obtained by the five algorithms on F5-F8. It is evident that the average MIGD values obtained by the five algorithms on the ZJZ problems are clearly worse than those obtained on the FDA and dMOP problems, in general. The high-severity change (nt = 5) greatly increases the difficulty of tracking the moving POS and POF for these algorithms, especially on F5. Compared with the competitors, MOEA/D-DM still performs best on most cases of these problems. MOEA/D-KF shows its superiority only on F5 with nt = 5, τt = 5 and 10. MOEA/D-DM has a slight advantage on F6 and F7 comparing with MOEA/D-KF. Note that although F8 is a triobjective optimization problem, its nonlinear correlation between decision variables is weaker than those of F5-F7. Consequently, the average MIGD values obtained by these algorithms on this problem are better than their results on F5-F7. And MOEA/D-DM still performs best among all algorithms, whereas MOEA/D and MOEA/D-KF perform similarly.

The tracking characteristics of these algorithms on ZJZ are presented in Fig. 6. We can see that the tracking curves of all algorithms on F5-F7 fluctuate widely in response to all environmental changes, due to the strong nonlinear characteristics of the varying POS and POF over time. The strong nonlinear correlation of POS between consecutive MOPs causes the prediction scheme difficulty in estimating relatively accurate new locations of the POS. Although the prediction scheme cannot estimate well on these problems, the three prediction-based algorithms (MOEA/DM, MOEA/D-KF and DDS) seem to have better tracking abilities on F5-F7 compared with the two partial random re-initialization-based algorithms (MOEA/D and DNSGA-II-A), in general. The weaker nonlinear correlation of the POS in F8 between consecutive MOPs decreases the difficulty of tracking the moving POS and POF. Consequently, the three prediction-based algorithms have much better tracking curve stability, and MOEA/D-DM is the best. The performances of these algorithms on F8 indicate that the prediction scheme works well on this problem, especially the difference model we proposed. 5.4. Results on JY The results obtained by the five algorithms on JY test suites are presented in Table 7. JY1 and JY2 are relatively simple problems, in which there is no correlation between decision variables and the

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

483

Table 4 Mean and standard deviations of MIGD values on FDA test suite obtained by five algorithms. Problems

(nt , τt )

MOEA/D-DM

MOEA/D

MOEA/D-KF

DNSGA-II-A

DDS

FDA1

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0128 ± 0.0007 0.0071 ± 0.0002 0.0049 ± 0.0001 0.0201 ± 0.0019 0.0082 ± 0.0002 0.0054 ± 0.0001

0.0297 ± 0.0007(+) 0.0116 ± 0.0002(+) 0.0065 ± 0.0001(+) 0.0807 ± 0.0122(+) 0.0169 ± 0.0007(+) 0.0072 ± 0.0001(+)

0.0249 ± 0.0024(+) 0.0099 ± 0.0001(+) 0.0060 ± 0.0001(+) 0.0392 ± 0.0161(+) 0.0101 ± 0.0002(+) 0.0060 ± 0.0001(+)

0.1135 ± 0.0044(+) 0.0381 ± 0.0020(+) 0.0173 ± 0.0005(+) 0.2005 ± 0.0053(+) 0.0792 ± 0.0031(+) 0.0249 ± 0.0010(+)

0.0596 ± 0.0041(+) 0.0236 ± 0.0011(+) 0.0104 ± 0.0002(+) 0.0757 ± 0.0036(+) 0.0392 ± 0.0021(+) 0.0161 ± 0.0011(+)

FDA2

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0089 ± 0.0006 0.0060 ± 0.0002 0.0048 ± 0.0001 0.0125 ± 0.0019 0.0070 ± 0.0007 0.0051 ± 0.0001

0.0107 ± 0.0005(+) 0.0070 ± 0.0004(+) 0.0052 ± 0.0001(+) 0.0157 ± 0.0011(+) 0.0084 ± 0.0005(+) 0.0057 ± 0.0001(+)

0.0587 ± 0.0179(+) 0.0107 ± 0.0006(+) 0.0063 ± 0.0002(+) 0.0573 ± 0.0204(+) 0.0110 ± 0.0004(+) 0.0061 ± 0.0001(+)

0.0221 ± 0.0016(+) 0.0119 ± 0.0007(+) 0.0087 ± 0.0001(+) 0.0342 ± 0.0014(+) 0.0154 ± 0.0006(+) 0.0099 ± 0.0004(+)

0.0312 ± 0.0041(+) 0.0126 ± 0.0004(+) 0.0105 ± 0.0012(+) 0.0541 ± 0.0039(+) 0.0190 ± 0.0012(+) 0.0103 ± 0.0003(+)

FDA3

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0552 ± 0.0074 0.0281 ± 0.0027 0.0155 ± 0.0007 0.0726 ± 0.0069 0.0471 ± 0.0067 0.0287 ± 0.0016

0.0736 ± 0.0088(+) 0.0511 ± 0.0084(+) 0.0302 ± 0.0068(+) 0.1089 ± 0.0070(+) 0.0671 ± 0.0060(+) 0.0510 ± 0.0049(+)

0.0732 ± 0.0099(+) 0.0466 ± 0.0073(+) 0.0370 ± 0.0086(+) 0.0795 ± 0.0097(+) 0.0533 ± 0.0006(+) 0.0349 ± 0.0033(+)

0.1256 ± 0.0061(+) 0.0960 ± 0.0055(+) 0.0812 ± 0.0017(+) 0.1668 ± 0.0072(+) 0.1016 ± 0.0054(+) 0.0902 ± 0.0046(+)

0.1007 ± 0.0126(+) 0.0852 ± 0.0029(+) 0.0781 ± 0.0091(+) 0.1031 ± 0.0109(+) 0.0882 ± 0.0071(+) 0.0852 ± 0.0019(+)

FDA4

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0469 ± 0.0005 0.0411 ± 0.0002 0.0380 ± 0.0001 0.0488 ± 0.0017 0.0411 ± 0.0006 0.0381 ± 0.0001

0.0572 ± 0.0007(+) 0.0454 ± 0.0002(+) 0.0397 ± 0.0006(+) 0.0576 ± 0.0011(+) 0.0439 ± 0.0001(+) 0.0393 ± 0.0002(+)

0.0544 ± 0.0013(+) 0.0430 ± 0.0004(+) 0.0391 ± 0.0003(+) 0.0545 ± 0.0004(+) 0.0429 ± 0.0001(+) 0.0388 ± 0.0005(+)

0.1257 ± 0.0023(+) 0.0715 ± 0.0005(+) 0.0543 ± 0.0007(+) 0.1816 ± 0.0047(+) 0.0949 ± 0.0010(+) 0.0569 ± 0.0007(+)

0.1376 ± 0.0080(+) 0.1269 ± 0.0095(+) 0.1262 ± 0.0119(+) 0.1527 ± 0.0027(+) 0.1501 ± 0.0088(+) 0.1288 ± 0.0102(+)

FDA5

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0986 ± 0.0073 0.0787 ± 0.0040 0.0611 ± 0.0042 0.1063 ± 0.0012 0.0818 ± 0.0032 0.0652 ± 0.0015

0.0982 ± 0.0050(−) 0.0705 ± 0.0069(−) 0.0670 ± 0.0056(+) 0.0992 ± 0.0016(−) 0.0656 ± 0.0013(−) 0.0551 ± 0.0015(−)

0.0854 ± 0.0153(−) 0.0602 ± 0.0009(−) 0.0568 ± 0.0063(−) 0.0873 ± 0.0012(−) 0.0636 ± 0.0073(−) 0.0574 ± 0.0050(−)

0.1049 ± 0.0010(+) 0.0608 ± 0.0070(−) 0.0539 ± 0.0031(−) 0.1387 ± 0.0013(+) 0.0767 ± 0.0067(−) 0.0471 ± 0.0042(−)

0.0813 ± 0.0017(−) 0.0725 ± 0.0015(−) 0.0680 ± 0.0050(+) 0.0841 ± 0.0085(−) 0.0819 ± 0.0072(∼) 0.0745 ± 0.0016(+)

Table 5 Mean and standard deviations of MIGD values on dMOP test suite obtained by five algorithms. Problems

(nt , τt )

MOEA/D-DM

MOEA/D

MOEA/D-KF

DNSGA-II-A

DDS

dMOP1

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0261 ± 0.0108 0.0081 ± 0.0030 0.0038 ± 0.0001 0.0298 ± 0.0100 0.0093 ± 0.0039 0.0039 ± 0.0001

0.0339 ± 0.0184(+) 0.0107 ± 0.0061(+) 0.0040 ± 0.0001(+) 0.0572 ± 0.0319(+) 0.0106 ± 0.0036(∼) 0.0047 ± 0.0009(+)

0.0303 ± 0.0102(+) 0.0125 ± 0.0025(+) 0.0042 ± 0.0001(+) 0.0594 ± 0.0297(+) 0.0104 ± 0.0070(∼) 0.0042 ± 0.0001(+)

0.1129 ± 0.0460(+) 0.0396 ± 0.0175(+) 0.0146 ± 0.0016(+) 0.1172 ± 0.0337(+) 0.0448 ± 0.0066(+) 0.0156 ± 0.0018(+)

0.0832 ± 0.0223(+) 0.0350 ± 0.0116(+) 0.0331 ± 0.0141(+) 0.0882 ± 0.0060(+) 0.0380 ± 0.0194(+) 0.0372 ± 0.0163(+)

dMOP2

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0261 ± 0.0087 0.0077 ± 0.0009 0.0047 ± 0.0001 0.0404 ± 0.0022 0.0129 ± 0.0024 0.0054 ± 0.0002

0.0449 ± 0.0014(+) 0.0141 ± 0.0026(+) 0.0063 ± 0.0001(+) 0.0914 ± 0.0075(+) 0.0203 ± 0.0006(+) 0.0069 ± 0.0001(+)

0.0319 ± 0.0105(+) 0.0098 ± 0.0004(+) 0.0055 ± 0.0001(+) 0.0416 ± 0.0062(+) 0.0135 ± 0.0005(+) 0.0057 ± 0.0001(+)

0.0982 ± 0.0026(+) 0.0407 ± 0.0005(+) 0.0140 ± 0.0006(+) 0.1637 ± 0.0068(+) 0.0656 ± 0.0015(+) 0.0198 ± 0.0005(+)

0.0538 ± 0.0060(+) 0.0256 ± 0.0081(+) 0.0195 ± 0.0072(+) 0.0638 ± 0.0126(+) 0.0431 ± 0.0036(+) 0.0390 ± 0.0161(+)

dMOP3

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0748 ± 0.0048 0.0184 ± 0.0006 0.0059 ± 0.0001 0.1138 ± 0.0091 0.0205 ± 0.0026 0.0064 ± 0.0001

0.0535 ± 0.0074(−) 0.0159 ± 0.0020(−) 0.0062 ± 0.0001(+) 0.0831 ± 0.0106(−) 0.0197 ± 0.0018(−) 0.0068 ± 0.0002(+)

0.0644 ± 0.0146(−) 0.0141 ± 0.0009(−) 0.0060 ± 0.0001(∼) 0.0845 ± 0.0131(−) 0.0156 ± 0.0091(−) 0.0062 ± 0.0001(∼)

0.0746 ± 0.0019(∼) 0.0277 ± 0.0011(+) 0.0107 ± 0.0002(+) 0.1193 ± 0.0035(+) 0.0478 ± 0.0016(+) 0.0147 ± 0.0003(+)

0.0809 ± 0.0055(+) 0.0443 ± 0.0030(+) 0.0307 ± 0.0104(+) 0.0913 ± 0.0107(−) 0.0639 ± 0.0037(+) 0.0490 ± 0.0163(+)

POS of JY1 and JY2 change over time in a sinusoidal pattern, despite the POF in JY2 changing its shape over time. These two problems are mainly used to test the convergence speed and the reactivity of an algorithm. As long as a DMOEA has fast convergence, it can easily solve these two problems. Accordingly, the three MOEA/Dbased algorithms perform better than the two NSGA-II-based algorithms on JY1 and JY2. Since the POS of JY1 and JY2 change over time in a weakly-nonlinear pattern, both the difference model proposed by us and the Kalman Filter proposed in MOEA/D-KF can improve the tracking ability of the algorithm. And our difference model seems better than the Kalman Filter model on these two problems, across different frequencies and severities of change. DDS performs better than DNSGA-II-A on JY1 and JY2, benefiting from their proposed prediction scheme that indeed improves the tracking ability of NSGA-II. JY3 introduces time-varying nonmonotonic dependencies between any two decision variables, meaning

that each variable has a different amount of change. This change does not obey the assumption that the solutions have the same motion as the centroid over time. Therefore, the difference model in our algorithm loses its superiority when dealing with this kind of problems. MOEA/D-KF shows its advantage on this problem, which may benefit from their use of Gaussian noise. But the superiority of MOEA/D-KF only appears when the environment changes slowly (τt = 10 and 20). JY4 is constructed to have a time-varying number of disconnected POF segments, (i.e. the POF is discontinuous). The difficulty of solving this problem for DMOEAs is that it is hard to cover all of the POF components, in addition to tracking the time-varying POF. As Table 7 shows, giving more evaluations in each time window cannot clearly improve these algorithms’ performances on JY4. The three MOEA/D-based algorithms actually perform similarly, and MOEA/D has a weak advantage on this problem. The POF of JY5 is simple and changes from convex

484

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Fig. 5. Tracking the IGD values with the environmental change obtained by five algorithms for FDA and dMOP problems with nt = 10, τt = 10, averaged over 30 runs.

geometry to concave geometry, whereas the POS of JY5 remains stationary. Accordingly, the strategy of retaining a portion of the old solutions obtained in the immediately previous environment shows its superiority. MOEA/D-DM and MOEA/D perform very well on this problem, even when the environment changes quickly and severely. JY6 is a multimodal problem, where the number of local optima changes over time, and the POS dynamically shifts in a sinusoidal pattern. Because of the existence of local optima,

MOEAs need more evaluations to approximate the POS well before the environment changes. Therefore, MOEA/D does not perform well until τt = 10 and 20, when nt = 10. The difference model in MOEA/D-DM still works better than the approach of partial random re-initialization in MOEA/D and the Kalman Filter model in MOEA/D-KF on this problem. JY7 is also a multimodal problem in which the number of local optima remains fixed, but the overall POF shape can be concave or convex due to environmental

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

485

Table 6 Mean and standard deviations of MIGD values on ZJZ test suite obtained by five algorithms. Problems

(nt , τt )

MOEA/D-DM

MOEA/D

MOEA/D-KF

DNSGA-II-A

DDS

F5

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.5565 ± 0.0514 0.1691 ± 0.0611 0.0518 ± 0.0094 1.0183 ± 0.2008 0.7266 ± 0.0596 0.5353 ± 0.0619

0.7451 ± 0.0695(+) 0.6148 ± 0.0597(+) 0.3610 ± 0.0672(+) 1.3587 ± 0.0694(+) 0.8078 ± 0.0864(+) 0.5622 ± 0.0525(+)

0.6644 ± 0.2153(+) 0.2738 ± 0.0690(+) 0.0891 ± 0.0431(+) 0.7054 ± 0.7055(−) 0.3495 ± 0.0596(−) 0.3918 ± 0.0605(−)

2.0327 ± 0.0764(+) 0.8145 ± 0.1088(+) 0.3482 ± 0.0058(+) 3.3076 ± 0.0957(+) 1.4451 ± 0.0763(+) 0.7117 ± 0.0442(+)

0.7313 ± 0.0285(+) 0.3515 ± 0.0213(+) 0.1463 ± 0.0037(+) 1.5468 ± 0.2493(+) 0.7022 ± 0.0788(−) 0.3189 ± 0.0157(−)

F6

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.3688 ± 0.0718 0.2304 ± 0.0550 0.1149 ± 0.0270 0.4232 ± 0.0618 0.4158 ± 0.0171 0.4085 ± 0.0147

0.6259 ± 0.1432(+) 0.2858 ± 0.0217(+) 0.1898 ± 0.0314(+) 0.8785 ± 0.1421(+) 0.6953 ± 0.0969(+) 0.4843 ± 0.0861(+)

0.4923 ± 0.1305(+) 0.2474 ± 0.0770(+) 0.1282 ± 0.0281(+) 0.7259 ± 0.0347(+) 0.5842 ± 0.0691(+) 0.4606 ± 0.0521(+)

1.4238 ± 0.0249(+) 0.7294 ± 0.0596(+) 0.2856 ± 0.0346(+) 2.1295 ± 0.0314(+) 1.1189 ± 0.0444(+) 0.6097 ± 0.0713(+)

1.0297 ± 0.1411(+) 0.3385 ± 0.0075(+) 0.1486 ± 0.0066(+) 1.8156 ± 0.3615(+) 0.8147 ± 0.0718(+) 0.4207 ± 0.0218(+)

F7

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.4115 ± 0.0338 0.1767 ± 0.0275 0.0997 ± 0.0221 0.5551 ± 0.0563 0.4418 ± 0.0196 0.3799 ± 0.0290

0.4735 ± 0.0529(+) 0.3590 ± 0.0352(+) 0.2544 ± 0.0400(+) 0.6335 ± 0.0559(+) 0.5227 ± 0.0751(+) 0.4086 ± 0.0331(+)

0.5940 ± 0.2637(+) 0.1993 ± 0.0207(+) 0.1223 ± 0.0473(+) 0.5889 ± 0.0968(+) 0.4971 ± 0.0249(+) 0.3971 ± 0.0371(+)

1.5068 ± 0.0382(+) 0.7448 ± 0.0615(+) 0.2842 ± 0.0173(+) 2.0911 ± 0.0531(+) 1.1363 ± 0.0689(+) 0.6095 ± 0.0988(+)

1.0893 ± 0.0843(+) 0.3621 ± 0.0335(+) 0.1479 ± 0.0111(+) 1.9872 ± 0.0827(+) 0.8015 ± 0.1077(+) 0.3992 ± 0.0192(+)

F8

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0869 ± 0.0056 0.0742 ± 0.0019 0.0569 ± 0.0024 0.1083 ± 0.0002 0.0868 ± 0.0028 0.0670 ± 0.0040

0.1212 ± 0.0032(+) 0.0966 ± 0.0052(+) 0.0717 ± 0.0008(+) 0.1423 ± 0.0019(+) 0.1061 ± 0.0038(+) 0.0839 ± 0.0007(+)

0.1432 ± 0.0009(+) 0.1081 ± 0.0045(+) 0.0787 ± 0.0052(+) 0.1478 ± 0.0005(+) 0.1029 ± 0.0041(+) 0.0841 ± 0.0006(+)

0.2648 ± 0.0031(+) 0.1568 ± 0.0037(+) 0.1301 ± 0.0002(+) 0.4124 ± 0.0016(+) 0.2231 ± 0.0025(+) 0.1364 ± 0.0057(+)

0.2847 ± 0.0008(+) 0.1820 ± 0.0050(+) 0.1439 ± 0.0112(+) 0.3066 ± 0.0199(+) 0.2267 ± 0.0006(+) 0.1756 ± 0.0140(+)

Fig. 6. Tracking curves obtained by the five algorithms on ZJZ test instances with nt = 10, τt = 10, averaged over 30 runs.

changes. As seen in Table 7, this problem is very challenging to all algorithms. The environment changes too fast for them to track the changed POF and POS. In JY8, the POS remains static, whereas the geometry and the number of mixed segments of the POF vary over time. These five algorithms seem to perform similarly on this problem, in general. MOEA/D-DM performs best when τt = 5, whereas DNSGA-II-A has better performance when τt = 10 or 20. The tracking characteristics of these algorithms on JY test instances are presented in Fig. 7. Since the POS of JY1 and JY2 change over time in a sinusoidal pattern, MOEA/D-DM can track the varying POS and POF effectively and stably on JY1 and JY2. MOEA/D

and MOEA/D-KF also show good tracking abilities on these two problems. Although the POS of JY3 changes over time in a stronglynonlinear pattern, MOEA/D-DM still has a stable tracking ability. The POS of JY4 also changes in a sinusoidal pattern; therefore, the three MOEA/D-based algorithms have stable tracking curves that also benefit from fast convergence, not only from the effect of a prediction scheme. Since the POS of JY5 does not change, MOEA/DDM and MOEA/D have the best and most stable tracking abilities, which may result from their retaining some old solutions in the new environment. The tracking curves of all algorithms fluctuate widely on JY6 and JY7, which illustrates that these two problems

486

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Table 7 Mean and standard deviations of MIGD values on JY test suite obtained by five algorithms. Problems

(nt , τt )

MOEA/D-DM

MOEA/D

MOEA/D-KF

DNSGA-II-A

DDS

JY1

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0119 ± 0.0005 0.0081 ± 0.0001 0.0070 ± 0.0001 0.0140 ± 0.0006 0.0088 ± 0.0001 0.0074 ± 0.0001

0.0221 ± 0.0013(+) 0.0112 ± 0.0004(+) 0.0079 ± 0.0001(+) 0.0331 ± 0.0024(+) 0.0124 ± 0.0002(+) 0.0080 ± 0.0001(+)

0.0165 ± 0.0007(+) 0.0095 ± 0.0001(+) 0.0074 ± 0.0001(+) 0.0168 ± 0.0008(+) 0.0094 ± 0.0002(+) 0.0075 ± 0.0001(∼)

2.2917 ± 0.1051(+) 2.2579 ± 0.2837(+) 2.1895 ± 0.0605(+) 2.3698 ± 0.1423(+) 2.3243 ± 0.0606(+) 2.2547 ± 0.0059(+)

0.1353 ± 0.0211(+) 0.0788 ± 0.0091(+) 0.0555 ± 0.0065(+) 0.1417 ± 0.0152(+) 0.0975 ± 0.0055(+) 0.0724 ± 0.0120(+)

JY2

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0119 ± 0.0004 0.0068 ± 0.0001 0.0054 ± 0.0001 0.0139 ± 0.0001 0.0075 ± 0.0001 0.0056 ± 0.0001

0.0211 ± 0.0001(+) 0.0101 ± 0.0002(+) 0.0062 ± 0.0001(+) 0.0321 ± 0.0019(+) 0.0114 ± 0.0004(+) 0.0065 ± 0.0001(+)

0.0153 ± 0.0006(+) 0.0082 ± 0.0002(+) 0.0058 ± 0.0001(+) 0.0165 ± 0.0006(+) 0.0082 ± 0.0001(+) 0.0059 ± 0.0001(+)

0.2006 ± 0.0077(+) 0.1114 ± 0.0104(+) 0.1216 ± 0.0132(+) 0.3141 ± 0.0028(+) 0.1809 ± 0.0119(+) 0.1274 ± 0.0071(+)

0.0934 ± 0.0085(+) 0.0514 ± 0.0019(+) 0.0289 ± 0.0009(+) 0.1197 ± 0.0123(+) 0.0761 ± 0.0080(+) 0.0414 ± 0.0050(+)

JY3

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.1014 ± 0.0002 0.0991 ± 0.0007 0.0977 ± 0.0002 0.1065 ± 0.0050 0.1001 ± 0.0002 0.0856 ± 0.0239

0.1026 ± 0.0009(+) 0.0987 ± 0.0003(−) 0.0977 ± 0.0001(∼) 0.1021 ± 0.0001(−) 0.1008 ± 0.0016(+) 0.0986 ± 0.0001(+)

0.1658 ± 0.0210(+) 0.0768 ± 0.0367(−) 0.0439 ± 0.0015(−) 0.1669 ± 0.0050(+) 0.0784 ± 0.0245(−) 0.0497 ± 0.0229(−)

0.1371 ± 0.0050(+) 0.1070 ± 0.0023(+) 0.1034 ± 0.0014(+) 0.1398 ± 0.0072(+) 0.1111 ± 0.0012(+) 0.1033 ± 0.0017(+)

0.1522 ± 0.0326(+) 0.1287 ± 0.0272(+) 0.0811 ± 0.0051(+) 0.1593 ± 0.0228(+) 0.1188 ± 0.0107(+) 0.1097 ± 0.0036(+)

JY4

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0703 ± 0.0001 0.0688 ± 0.0002 0.0646 ± 0.0004 0.0697 ± 0.0002 0.0677 ± 0.0001 0.0634 ± 0.0003

0.0693 ± 0.0001(−) 0.0662±0.0002(−) 0.0616 ± 0.0009(−) 0.0689 ± 0.0001(−) 0.0648 ± 0.0001(−) 0.0633 ± 0.0011(∼)

0.0690 ± 0.0003(−) 0.0665 ± 0.0002(−) 0.0623 ± 0.0004(−) 0.0690 ± 0.0001(−) 0.0665 ± 0.0005(−) 0.0637 ± 0.0018(+)

0.3046 ± 0.0350(+) 0.2134 ± 0.0078(+) 0.2018 ± 0.0181(+) 0.4135 ± 0.0148(+) 0.2791 ± 0.0182(+) 0.2725 ± 0.0238(+)

0.1143 ± 0.0022(+) 0.0834 ± 0.0037(+) 0.0635 ± 0.0016(−) 0.1322 ± 0.0059(+) 0.1164 ± 0.0126(+) 0.0792 ± 0.0049(+)

JY5

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0082 ± 0.0002 0.0070 ± 0.0002 0.0067 ± 0.0001 0.0087 ± 0.0001 0.0072 ± 0.0001 0.0068 ± 0.0001

0.0086 ± 0.0003(+) 0.0072 ± 0.0001(+) 0.0069 ± 0.0001(+) 0.0088 ± 0.0001(∼) 0.0073 ± 0.0001(+) 0.0069 ± 0.0001(∼)

0.0185 ± 0.0003(+) 0.0099 ± 0.0002(+) 0.0076 ± 0.0001(+) 0.0189 ± 0.0004(+) 0.0101 ± 0.0002(+) 0.0076 ± 0.0001(+)

0.2129 ± 0.0163(+) 0.2098 ± 0.0408(+) 0.2063 ± 0.0079(+) 0.1037 ± 0.0265(+) 0.0866 ± 0.0073(+) 0.0670 ± 0.0079(+)

0.0325 ± 0.0030(+) 0.0183 ± 0.0021(+) 0.0132 ± 0.0002(+) 0.0355 ± 0.0051(+) 0.0199 ± 0.0015(+) 0.0140 ± 0.0021(+)

JY6

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.7683 ± 0.2305 0.0709 ± 0.0251 0.0135 ± 0.0003 1.7921 ± 0.0977 1.0961 ± 0.0252 0.3247 ± 0.0559

1.2888 ± 0.1813(+) 0.5860 ± 0.0408(+) 0.2530 ± 0.0079(+) 2.0366 ± 0.0671(+) 1.1478 ± 0.1417(+) 0.4961 ± 0.0239(+)

2.0337 ± 0.2023(+) 0.4830 ± 0.1970(+) 0.1443 ± 0.0197(+) 2.5656 ± 0.1670(+) 1.2035 ± 0.0854(+) 0.3361 ± 0.0701(+)

2.9309 ± 0.1072(+) 1.8941 ± 0.0671(+) 1.1829 ± 0.0125(+) 3.6841 ± 0.0481(+) 2.4085 ± 0.0481(+) 1.6585 ± 0.0550(+)

2.0209 ± 0.0813(+) 1.5976 ± 0.0462(+) 1.3639 ± 0.0307(+) 2.6077 ± 0.1800(+) 2.1725 ± 0.0874(+) 2.0815 ± 0.0340(+)

JY7

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

2.0947 ± 0.2367 1.5058 ± 0.3002 1.0383 ± 0.1342 2.0939 ± 0.3763 1.9093 ± 0.2225 1.1910 ± 0.1973

2.1384 ± 0.2570(+) 2.0049 ± 0.5105(+) 1.5039 ± 0.6146(+) 2.2279 ± 0.2962(+) 1.2515 ± 0.1999(−) 1.1790 ± 0.2316(−)

3.1103 ± 0.9544(+) 3.8106 ± 0.6338(+) 2.4452 ± 0.4482(+) 2.6165 ± 0.3025(+) 2.0159 ± 0.6514(+) 1.6991 ± 0.7482(+)

8.3470 ± 0.2555(+) 5.0254 ± 0.4324(+) 1.6601 ± 0.1095(+) 9.2693 ± 0.4685(+) 6.9805 ± 0.3010(+) 3.6907 ± 0.1779(+)

5.3513 ± 0.3697(+) 3.8062 ± 0.4774(+) 2.6459 ± 0.1799(+) 5.6526 ± 0.6497(+) 3.7794 ± 0.3253(+) 2.3129 ± 0.1242(+)

JY8

(10,5) (10,10) (10,20) (5,5) (5,10) (5,20)

0.0226 ± 0.0012 0.0189 ± 0.0007 0.0186 ± 0.0002 0.0256 ± 0.0012 0.0209 ± 0.0006 0.0204 ± 0.0004

0.0236 ± 0.0011(+) 0.0204 ± 0.0012(+) 0.0195 ± 0.0008(+) 0.0258 ± 0.0008(+) 0.0220 ± 0.0009(+) 0.0219 ± 0.0024(+)

0.0330 ± 0.0027(+) 0.0219 ± 0.0008(+) 0.0213 ± 0.0025(+) 0.0343 ± 0.0008(+) 0.0249 ± 0.0017(+) 0.0242 ± 0.0007(+)

0.0330 ± 0.0039(+) 0.0153 ± 0.0012(−) 0.0121 ± 0.0012(−) 0.0370 ± 0.0036(+) 0.0167 ± 0.0015(−) 0.0124 ± 0.0014(−)

0.0706 ± 0.0016(+) 0.0257 ± 0.0030(+) 0.0156 ± 0.0012(−) 0.1734 ± 0.0158(+) 0.0508 ± 0.0066(+) 0.0189 ± 0.0012(−)

greatly challenge DMOEAs because of the complex variation of their POF. DNSGA-II-A has relatively stable tracking ability on JY8 compared with the other algorithms, especially after the fifth time window. 5.5. Further analysis The above experimental studies have clearly showed the general, although not universal, superiority of the proposed MOEA/DDM for handling DMOPs. This section will analyze some key components in the algorithm. (1) Parameter sensitivity: In MOEA/D-DM, we use a parameter q (q ∈ [0, 1]) to control the number of predicted solutions in the new population when an environmental change is detected. A smaller q means that more old solutions obtained in the immediately previous environment will be retained in the new population, whereas a larger q means more predicted solutions will be placed in the new population. A proper value of q is determined by the shift distance of the POS in the decision space between two sequential

time windows. If the location of the POS in the later time window is very close to the location in the immediately previous time window, more old solutions should be retained. Otherwise, more retained old solutions may not be useful in the new environment. The default value of q in this paper is 0.5, and now we will explore how varying this value affects the performance of MOEA/D-DM. We choose FDA1, JY2, F5–F6, dMOP1 and JY5 to investigate the influence of varying the value of q. The POS and/or POF of FDA1 and JY2 change over time in weakly-nonlinear patterns, whereas the POS and POF of F5–F6 change over time in strongly-nonlinear patterns. The POS of dMOP1 and JY5 remain static during all environmental changes. Three different severities of environmental change (nt = 5, 10 and 20) are used here, with τt = 10. Mean and standard deviations of MIGD values obtained by MOEA/D-DM with different values of q are presented in Table 8, in which the best value in each test instance is bolded. Based on the results on FDA1 and JY2, retaining all old solutions (q = 0) cannot obtain good results compared with the performance of MOEA/DDM with other values of q. Placing predicted solutions in the new population is obviously effective. The performance of MOEA/D-DM

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

487

Fig. 7. Tracking curves obtained by the five algorithms on JY test instances with nt = 10, τt = 10, averaged over 30 runs.

is not sensitive to the value of q (as long as q ̸ = 0) on these two problems in which the POS and/or POF change over time in weaklynonlinear patterns. Similar to solving FDA1 and JY2, retaining only old solutions (q = 0) also cannot obtain good results when solving F5 and F6. Placing all predicted solutions in the population (q = 1) can obtain the best result when solving F5 with nt = 5; however the superiority is not retained when the severity of change becomes slight. Placing all predicted solutions in the population (q = 1) performs even worse than the approach of placing some

predicted solutions in the population (q = 0.3, 0.5 and 0.7) when solving F6 with nt = 5 and 10. When solving dMOP1, in which the POS does not change, retaining all old solutions or accepting all predicted solutions both seem improper if nt = 5 or 10, due to the large shift of the POF in the objective space. The mixture of old solutions and predicted solutions is more proper. When the POF changes very slightly (nt = 20), retaining only old solutions performs best. MOEA/D-DM with different values of q performs similarly on JY5 with different severities of changes. The results of

488

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

Table 8 Mean and standard deviations of MIGD values obtained by MOEA/D-DM with different values of q. Problems

(nt , τt )

q=0

q = 0.3

q = 0.5

q = 0.7

q=1

FDA1

(5,10) (10,10) (20,10)

0.0183 ± 0.0021 0.0131 ± 0.0031 0.0091 ± 0.0001

0.0088 ± 0.0002 0.0077 ± 0.0004 0.0070 ± 0.0001

0.0082 ± 0.0002 0.0071 ± 0.0002 0.0064 ± 0.0001

0.0079 ± 0.0001 0.0068 ± 0.0001 0.0067 ± 0.0003

0.0079 ± 0.0003 0.0071 ± 0.0005 0.0071 ± 0.0009

JY2

(5,10) (10,10) (20,10)

0.0113 ± 0.0003 0.0101 ± 0.0002 0.0088 ± 0.0001

0.0078 ± 0.0001 0.0073 ± 0.0001 0.0069 ± 0.0001

0.0075 ± 0.0001 0.0068 ± 0.0001 0.0065 ± 0.0001

0.0076 ± 0.0001 0.0070 ± 0.0001 0.0064 ± 0.0001

0.0073 ± 0.0001 0.0065 ± 0.0001 0.0061 ± 0.0001

F5

(5,10) (10,10) (20,10)

0.9736 ± 0.0632 0.8185 ± 0.1263 0.3543 ± 0.0824

0.7781 ± 0.1417 0.1861 ± 0.0837 0.0819 ± 0.0333

0.7266 ± 0.0596 0.1691 ± 0.0612 0.0716 ± 0.0177

0.6678 ± 0.1360 0.2116 ± 0.1218 0.0539 ± 0.0040

0.4341 ± 0.1260 0.2983 ± 0.1329 0.0689 ± 0.0153

F6

(5,10) (10,10) (20,10)

0.7625 ± 0.0027 0.2876 ± 0.0342 0.2608 ± 0.0552

0.4338 ± 0.0732 0.2588 ± 0.0691 0.0620 ± 0.0005

0.4158 ± 0.0171 0.2304 ± 0.0550 0.0539 ± 0.0059

0.4110 ± 0.0036 0.2716 ± 0.0171 0.0525 ± 0.0057

0.6379 ± 0.0088 0.4291 ± 0.0007 0.0480 ± 0.0095

dMOP1

(5,10) (10,10) (20,10)

0.0145 ± 0.0063 0.0129 ± 0.0001 0.0044 ± 0.0003

0.0082 ± 0.0070 0.0071 ± 0.0011 0.0068 ± 0.0031

0.0093 ± 0.0039 0.0081 ± 0.0030 0.0070 ± 0.0052

0.0095 ± 0.0058 0.0086 ± 0.0012 0.0079 ± 0.0003

0.0101 ± 0.0052 0.0131 ± 0.0048 0.0085 ± 0.0005

JY5

(5,10) (10,10) (20,10)

0.0075 ± 0.0001 0.0071 ± 0.0001 0.0072 ± 0.0002

0.0072 ± 0.0001 0.0070 ± 0.0002 0.0071 ± 0.0001

0.0072 ± 0.0001 0.0070 ± 0.0001 0.0071 ± 0.0001

0.0074 ± 0.0001 0.0071 ± 0.0001 0.0073 ± 0.0001

0.0075 ± 0.0003 0.0076 ± 0.0001 0.0077 ± 0.0001

Table 8 illustrate that the mixture of old solutions and predicted solutions is generally better than either approach alone, and the performance of MOEA/D-DM is robust across a variety of value of q (q ∈ (0, 1)). (2) Effect of nr in Replacement: In the original MOEA/D algorithm, the maximum number of solutions replaced by a child solution is bounded by nr , which was recommended to be set to be much smaller than the neighborhood size, in order to maintain the population diversity [32]. However, solving DMOPs requires that DMOEAs have fast convergence not only population diversity. A larger nr allows a high-quality child solution to replace many current solutions in its neighboring subproblems, which is called global replacement [49]. As a result, convergence is promoted, whereas diversity is decreased. This section investigates the effect of nr on the performance of MOEA/D-DM. In the above experiments, the default value of nr is set to be |E |. Now we use the recommended value nr = 2 in [32] to investigate the performance of MOEA/D-DM. FDA1, FDA4, dMOP1, F5, JY3 and JY6 are chosen to investigate the effect of two values of nr on MOEA/D-DM. There is no correlation between decision variables of FDA1, FDA4 or dMOP1; convergence is more needed than diversity for these problems. However, there is a strongly-nonlinear correlation among the decision variables of F5, which requires a diverse population during searching in order to find the entire POF. JY3 introduces timevarying nonmonotonic dependencies between any two decision variables, which tests the diversity performance of an algorithm in a dynamic environment. JY6 is a multimodal problem, which also requires diversity of the population to find the entire POF. Three different frequencies of environmental change (τt = 5, 10 and 20) are used, with nt = 10. Mean and standard deviations of MIGD values obtained by MOEA/D-DM with different values of nr are presented in Table 9, in which the best value in each test instance is bolded. Obviously, MOEA/D-DM with local replacement (nr = 2) performs worse than with global replacement (nr = |E |) on each kind of DMOP, especially on problems that change frequently. Although solving F5, JY3 and JY6 requires the population to maintain good diversity during optimization in a specific environment, the global replacement can balance the diversity and convergence to obtain better results. The results of Table 9 illustrate that the global replacement is more proper for solving DMOPs, due to its good tradeoff between convergence and diversity. (3) Study of different difference models in MOEA/D-DM: In our MOEA/D-DM, we use the second-order difference model starting

from the fourth time window to predict the new locations of the POS in the new environment, whereas the first-order difference model is applied in the third time window. This section is devoted to studying the effects of using different difference models in MOEA/D-DM. First we test using only the first-order difference model (denoted as MOEA/D-FDM), starting from the third time window, to predict the new locations of the POS in the new environment. Theoretically, the first-order difference model is suitable for those problems in which the POS changes over time in linear patterns or weakly-nonlinear patterns with slight severity. We choose FDA1, JY1, JY2, F5–F6 and JY3 to investigate performances of these two difference models. The POS and/or POF of FDA1, JY1 and JY2 change over time in weakly-nonlinear patterns, whereas the POS and POF of F5–F6 and JY3 change over time in stronglynonlinear patterns. The environment changes with different severities (nt = 3, 5 and 10) and a fixed frequency (τt = 10). Mean and standard deviations of MIGD values obtained by MOEA/D-DM and MOEA/D-FDM are presented in Table 10, in which the best value in each test instance is bolded. Results on FDA1, JY1 and JY2 show that the second-order difference model used in the algorithm is better than the first-order difference model used when the problems change severely (nt = 3 and 5). Although these problems change in weakly-nonlinear patterns, the POS still shifts a lot in the decision space when the environment changes severely. The new locations of the POS of these problems predicted by the first-order difference model may be further away from the true locations than the locations predicted by the second-order difference model. When the environments in these problems change slightly (nt = 10), these two models perform similarly. Results of F5 and F6 show that the second-order difference model performs much better than the first-order difference model, especially when the problems change severely. The two models perform similarly on JY3, since the dependencies between any two decision variables of JY3 become increasingly complicated over time, making it hard for the prediction model to estimate the new locations of the POS. Neither the second-order nor the first-order model can estimate well, although the environment changes only slightly. The results of Table 10 illustrate that the second-order difference model is superior when solving the problems that change severely in weaklynonlinear patterns or change in strongly-nonlinear patterns. 6. Discussion and conclusions This paper presents a simple yet effective evolutionary dynamic multiobjective optimization algorithm based on decomposition.

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

489

Table 9 Mean and standard deviations of MIGD values obtained by MOEA/D-DM with different values of nr . Problems

(nt , τt )

nr = |E |

nr = 2

FDA1

(10,5) (10,10) (10,20)

0.0128 ± 0.0007 0.0071 ± 0.0002 0.0049 ± 0.0001

0.0513 ± 0.0006 0.0295 ± 0.0009 0.0176 ± 0.0007

FDA4

(10,5) (10,10) (10,20)

0.0469 ± 0.0005 0.0411 ± 0.0002 0.0380 ± 0.0001

0.1042 ± 0.0003 0.0707 ± 0.0005 0.0592 ± 0.0014

dMOP1

(10,5) (10,10) (10,20)

0.0261 ± 0.0108 0.0081 ± 0.0030 0.0039 ± 0.0001

0.0272 ± 0.0040 0.0112 ± 0.0012 0.0066 ± 0.0004

F5

(10,5) (10,10) (10,20)

0.5565 ± 0.0514 0.1691 ± 0.0612 0.0518 ± 0.0094

0.6425 ± 0.0801 0.2857 ± 0.1172 0.0864 ± 0.0517

JY3

(10,5) (10,10) (10,20)

0.1014 ± 0.0002 0.0991 ± 0.0007 0.0977 ± 0.0002

0.1316 ± 0.0023 0.1121 ± 0.0023 0.1046 ± 0.0002

JY6

(10,5) (10,10) (10,20)

0.7683 ± 0.2305 0.0709 ± 0.0251 0.0135 ± 0.0003

2.1457 ± 0.0837 1.3849 ± 0.0378 0.6151 ± 0.1080

Table 10 Mean and standard deviations of MIGD values obtained by MOEA/D-DM with different models. Problems

(nt , τt )

MOEA/D-DM

MOEA/D-FDM

FDA1

(3,10) (5,10) (10,10)

0.0107 ± 0.0003 0.0082 ± 0.0002 0.0071 ± 0.0002

0.0131 ± 0.0009 0.0096 ± 0.0003 0.0072 ± 0.0002

JY1

(3,10) (5,10) (10,10)

0.0102 ± 0.0005 0.0088 ± 0.0001 0.0081 ± 0.0001

0.0110 ± 0.0002 0.0099 ± 0.0001 0.0086 ± 0.0001

JY2

(3,10) (5,10) (10,10)

0.0092 ± 0.0002 0.0075 ± 0.0001 0.0068 ± 0.0001

0.0097 ± 0.0001 0.0085 ± 0.0001 0.0069 ± 0.0001

F5

(3,10) (5,10) (10,10)

1.0445 ± 0.1083 0.7266 ± 0.0596 0.1691 ± 0.0612

1.5111 ± 0.0614 0.7968 ± 0.0157 0.4645 ± 0.2579

F6

(3,10) (5,10) (10,10)

0.7252 ± 0.1137 0.4158 ± 0.0171 0.2304 ± 0.0550

0.9050 ± 0.0627 0.7750 ± 0.0316 0.2935 ± 0.0466

JY3

(3,10) (5,10) (10,10)

0.1021 ± 0.0004 0.1001 ± 0.0002 0.0991 ± 0.0007

0.1031 ± 0.0002 0.1003 ± 0.0002 0.0991 ± 0.0003

The algorithm uses a difference model to predict the new locations of the POS when detecting an environmental change. Unlike many prediction schemes in the literature, we did not predict each solution’s new location directly. In our difference model, the motion of the entire POS is represented by the trajectory of the centroid. The model uses the two or three most recent previous locations of the centroid to estimate its later motion, including the direction and step size of the motion. The other members of the POS are assumed to have the same motion as the centroid at any time, and their new locations are predicted based on their current locations and the estimated motion. The predicted solutions are combined with a subset of the old solutions to explore in the new environment. MOEA/D-DM was compared with four state-of-the-art dynamic MOEAs through different kinds of benchmark functions. The experimental results indicated that our proposed MOEA/D-DM clearly outperforms the other algorithms on the majority of the benchmark functions, except on some problems that change over time in strongly-nonlinear patterns and with other characteristics that we have analyzed in the above experimental studies. The main novelty and advantage of this algorithm is that the first-order difference model (also called the linear prediction model in some literature) has been extended to a second-order difference model. Therefore, more historical solutions can be used to estimate the new locations in the new environment. Note that, although

a second-order linear model is used in MOEA/D-KF, the secondorder vector is controlled by the Gaussian noise and the previous locations are not involved in updating the new second-order vector, as can be deduced from the state-transition model. Another advantage of MOEA/D-DM is that there is only one additional parameter (q) added to the original MOEA/D, and the algorithm is not very sensitive to this parameter when solving DMOPs. In addition, we do not need to update this parameter over time, so it is much easier to implement this prediction model comparing with a Kalman filter. The accuracy of predicting solutions in our proposed model is more precise on most benchmark problems compared to the other algorithms, which can be seen in the experimental results. The excellent tracking ability of MOEA/D-DM also benefits from the global replacement used, which greatly promotes convergence speed of the algorithm and keeps a good tradeoff between population diversity and convergence. This proposed algorithm focuses on those problems that change intermittently over time according to relatively low-order rules, and shows great superiority on solving these problems. The proposed difference model is based on an assumption that each Pareto-optimal solution has the same motion as the centroid over time, and the changes in all tested problems except JY3 obey this assumption. More problems that go against this assumption should be tested to see whether this strategy still works. In addition, the

490

L. Cao, L. Xu, E.D. Goodman et al. / Applied Soft Computing Journal 76 (2019) 473–490

performance of MOEA/D-DM on a broader class of problems should be investigated—e.g., problems changing over time according to other nonlinear functions, not limited to the trigonometric functions used in some problems explored here. It may be that decision variables of problems changing according to different types of rules may, in fact, be closer to many real-world dynamic optimization problems. Acknowledgments This work was supported in part by the National Natural Science Foundation of China under Grant 61573258, and in part by U.S. National Science Foundation’s BEACON Center for the Study of Evolution in Action, funded under Cooperative Agreement No. DBI0939454. References [1] R. Azzouz, S. Bechikh, L. Ben Said, Multi-objective optimization with dynamic constraints and objectives : new challenges for evolutionary algorithms, in: Proceedings of ACM GECCO Conference, 2015, pp. 615–622. [2] A. Zhou, B.-Y. Qu, H. Li, S.-Z. Zhao, P.N. Suganthan, Q. Zhang, Multiobjective evolutionary algorithms: A survey of the state of the art, Swarm Evol. Comput. 1 (1) (2011) 32–49. [3] R. Azzouz, S. Bechikh, L. Ben Said, A multiple reference point-based evolutionary algorithm for dynamic multi-objective optimization with undetectable changes, in: Proc.IEEE CEC, 2014, pp. 3168–3175. [4] M. Helbig, A.P. Engelbrecht, Population-based metaheuristics for continuous boundary-constrained dynamic multi-objective optimisation problems, Swarm Evol. Comput. 14 (2014) 31–47. [5] R. Azzouz, S. Bechikh, L. Ben Said, Dynamic multi-objective optimization using evolutionary algorithms: a survey, in: S. Bechikh, R. Datta, A. Gupta (Eds.), Recent Adavances in Evolutionary Multi-objective Opimitzation.Adaptation, Learning, and Optimization, Vol. 20, Springer International Publishing Switzerland, 2017, pp. 31–70. [6] G.G. Yen, H. Lu, Dynamic multiobjective evolutionary algorithm: adaptive cell-based rank and density estimation, IEEE Trans. Evol. Comput. 7 (3) (2003) 253–274. [7] Y. Jin, J. Branke, Evolutionary optimization in uncertain environments-a survey, IEEE Trans. Evol. Comput. 9 (3) (2005) 303–317. [8] E. Tantar, A.-A. Tantar, P. Bouvry, On dynamic multi-objective optimization, classification and performance measures, in: Proc.IEEE CEC, 2011, pp. 2759– 2766. [9] R. Chen, K. Li, X. Yao, Dynamic multi-objectives optimization with a changing number of objectives, IEEE Trans. Evol. Comput. XX (X) (2017) 1–15. [10] S.U. Guan, Q. Chen, W. Mo, Evolving dynamic multi-objective optimization problems with objective replacement, Artif. Intell. Rev. 23 (3) (2005) 267–293. [11] A. Muruganantham, K.C. Tan, P. Vadakkepat, Evolutionary dynamic multiobjective optimization via kalman filter prediction, IEEE Trans. Cybern. 46 (12) (2016) 2862–2873. [12] S. Jiang, S. Yang, A steady-state and generational evolutionary algorithm for dynamic multiobjective optimization, IEEE Trans. Evol. Comput. 21 (1) (2017) 65–82. [13] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput. 6 (2) (2002) 182–197. [14] Q. Zhang, H. Li, MOEA/D: a multiobjective evolutionary algorithm based on decomposition, IEEE Trans. Evol. Comput. 11 (6) (2007) 712–731. [15] W.K. Mashwani, A. Salhi, M.A. Jan, M. Sulaiman, R.A. Khanum, A. Algarni, Evolutionary algorithms based on decomposition and indicator functions: state-of-the-art survey, Int. J. Adv. Comput. Sci. Appl. 7 (2) (2016) 583–593. [16] Q. Zhang, A. Zhou, Y. Jin, RM-MEDA: a regularity model based multiobjective estimation of distribution algorithm, IEEE Trans. Evol. Comput. 12 (1) (2008) 41–63. [17] C.K. Goh, K.C. Tan, A competitive-cooperative coevolutionary paradigm for dynamic multiobjective optimization, IEEE Trans. Evol. Comput. 13 (1) (2009) 103–127. [18] M. Helbig, A.P. Engelbrecht, Analysing the performance of dynamic multiobjective optimisation algorithms, in: Proc.IEEE CEC, 2013, pp. 1531–1539. [19] R. Liu, J. Fan, L. Jiao, Integration of improved predictive model and adaptive differential evolution based dynamic multi-objective evolutionary optimization algorithm, Appl. Intell. (2015) 192–207. [20] J. Wei, L. Jia, A novel particle swarm optimization algorithm with local search for dynamic constrained multi-objective optimization problems, in: Proc.IEEE CEC, 2013, pp. 2436–2443. [21] B.K. Deb, N. Rao, S. Karthik, Dynamic multi-objective optimization and decision-making using modified NSGA-II: a case study on hydro-thermal power scheduling, in: Proceedings of EMO, in: LNCS 4403, 2007, pp. 803–817.

[22] Z. Bingul, Adaptive genetic algorithms applied to dynamic multiobjective problems, Appl. Soft Comput. J. 7 (3) (2007) 791–799. [23] S. Sahmoud, H.R. Topcuoglu, a Memory-based NSGA-II; Algorithm for Dynamic Multi-objective Optimization Problems, in: EvoApplications, in: LNCS 9598, vol. 9598, 2016, pp. 296–310. [24] X. Chen, D. Zhang, X. Zeng, A stable matching-based selection and memory enhanced MOEA/D for evolutionary dynamic multiobjective optimization, in: Proceedings of IEEE 27th International Conference on Tools with Artificial Intelligence, 2015, pp. 478–485. [25] Y. Wang, B. Li, Investigation of memory-based multi-objective optimization evolutionary algorithm in dynamic environment, in: Proceeding of IEEE Congress on Evolutionary Computation, 2009, pp. 630–637. [26] A. Zhou, Y. Jin, Q. Zhang, B. Sendhoff, E. Tsang, Prediction-Based population re-initialization for evolutionary dynamic multi-objective optimization, in: Proceedings of EMO, in: LNCS4403, 2007, pp. 832–846. [27] A. Zhou, Y. Jin, Q. Zhang, A population prediction strategy for evolutionary dynamic multiobjective optimization, IEEE Trans. Cybern. 44 (1) (2014) 40– 53. [28] I. Hatzakis, D. Wallace, Dynamic multi-objective optimization evolutionary algorithms: a Forward-Looking approach, in: Proceedings of GECCO, Vol. 4, 2006, pp. 1201–1208. [29] I. Hatzakis, D. Wallace, Topology of anticipatory populations for evolutionary dynamic multi-objective optimization, in: 11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2006, pp. 1–10. [30] W.T. Koo, C.K. Goh, K.C. Tan, A predictive gradient strategy for multiobjective evolutionary algorithms in a fast changing environment, Memetic Comput. 2 (2) (2010) 87–110. [31] M. Farina, K. Deb, P. Amato, Dynamic multiobjective optimization problem: test cases, approximation, and applications, IEEE Trans. Evol. Comput. 8 (5) (2004) 425–442. [32] H. Li, Q. Zhang, Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II, IEEE Trans. Evol. Comput. 13 (2) (2009) 284–302. [33] W. Luo, J. Sun, C. Bu, H. Liang, Species-based particle swarm optimizer enhanced by memory for dynamic optimization, Appl. Soft Comput. 47 (2016) 130–140. [34] A.M. Turky, S. Abdullah, A multi-population harmony search algorithm with external archive for dynamic optimization problems, Inf. Sci. (Ny). 272 (2014) 84–95. [35] R. Liu, Y. Chen, W. Ma, C. Mu, L. Jiao, A novel cooperative coevolutionary dynamic multi-objective optimization algorithm using a new predictive model, Soft Comput. 18 (10) (2014) 1913–1929. [36] Y. Wu, Y. Jin, X. Liu, A directed search strategy for evolutionary dynamic multiobjective optimization, Soft Comput. 19 (11) (2015) 3221–3235. [37] A. Muruganantham, Y. Zhao, S.B. Gee, X. Qiu, K.C. Tan, Dynamic multiobjective optimization using evolutionary algorithm with kalman filter, Procedia Comput. Sci. 24 (2013) 66–75. [38] K.R. Harrison, B.M. Ombuki-Berman, A.P. Engelbrecht, Dynamic multiobjective optimization using charged vector evaluated particle swarm optimization, in: Proc.IEEE CEC, 2014, pp. 1929–1936. [39] M. Greeff, A.P. Engelbrecht, Solving dynamic multi-objective problems with vector evaluated particle swarm optimisation, in: Proc.IEEE CEC, 2008, pp. 2917–2924. [40] W.K. Mashwani, A. Salhi, A decomposition-based hybrid multiobjective evolutionary algorithm with dynamic resource allocation, Appl. Soft Comput. 12 (9) (2012) 2765–2780. [41] W.K. Mashwani, A. Salhi, O. Yeniay, M.A. Jan, R.A. Khanum, Hybrid adaptive evolutionary algorithm based on decomposition, Appl. Soft Comput. 57 (2017) 363–378. [42] W.K. Mashwani, A. Salhi, Multiobjective memetic algorithm based on decomposition, Appl. Soft Comput. 21 (2014) 221–243. [43] S. Jiang, S. Yang, Evolutionary dynamic multiobjective optimization: benchmarks and algorithm comparisons, IEEE Trans. Cybern. 47 (99) (2017) 198– 211. [44] W.K. Mashwani, A. Salhi, Multiobjective evolutionary algorithm based on multimethod with dynamic resources allocation, Appl. Soft Comput. 39 (2016) 292–309. [45] W.K. Mashwani, A. Salhi, O. Yeniay, H. Hussian, M.A. Jan, Hybrid nondominated sorting genetic algorithm with adaptive operators selection, Appl. Soft Comput. 56 (2017) 1–18. [46] M. Helbig, A.P. Engelbrecht, Performance measures for dynamic multiobjective optimisation algorithms, Inf. Sci. (Ny). 250 (2013) 61–81. [47] C.K. Goh, K.C. Tan, An investigation on noisy environments in evolutionary multiobjective optimization, IEEE Trans. Evol. Comput. 11 (3) (2007) 354–381. [48] X. Li, J. Branke, M. Kirley, On performance metrics and particle swarm methods for dynamic multiobjective optimization problems, in: Proc.IEEE CEC, 2007, pp. 576–583. [49] Z. Wang, Q. Zhang, A. Zhou, M. Gong, L. Jiao, Adaptive replacement strategies for MOEA/D, IEEE Trans. Cybern. 46 (2) (2016) 474–486.