A dynamic multi-objective evolutionary algorithm using center and multi-direction prediction strategies

1.   Introduction

A large number of multi-objective optimizations which change over time can be found in scientific and engineering fields, which are always known as dynamic multi-objective optimization problems (DMOPs) ^[1]. When the state changes over time, some multi-objective optimization models in path planning ^[2], resource distribution ^[3], energy scheduling ^[4,5], network layout ^[6], and electrode-magnetic micro-mirrors become dynamic ones ^[7]. Therefore, it is necessary to study the optimal solutions to DMOPs at each moment as well as convergence. More and more researchers devote themselves to solving such problems in recent years. In any current environment of dynamic optimization, the Pareto optimal solutions (PS) and fronts (PF) need to be found in time before the next environment is detected, which brings great challenges for optimization approaches to deal with this kind of problem.

The types of DMOPs can be found in ^[8,9]. A general mode of DMOPs can be formulated as follows:

where $m$ is the scale of objectives. The objective $f\left(x, t\right)$ is a function on the decision-making vector, which can change over time or environment $t$ taken from the set $\left\{0, \mathrm{ }1, \mathrm{ }2, \cdots, T\right\}$. $T$ is the total number of environment changes, and the objective space is denoted by ${R}^{m}.$ Decision-making vector $x = \left({x}_{1}\left(t\right), {x}_{2}\left(t\right), \cdots, {x}_{n}\left(t\right)\right)$ is n-dimensional in space $\Omega.$ Both the decision variable ${x}_{i}\left(t\right), i = 1, 2, ..., n,$ and objective function ${f}_{k}\left(x, t\right), k = 1, 2, ..., m,$ depend on time variable $t$ and dynamically change, once the new environment is detected.

DMOPs consist of static multi-objective optimization problems (MOPs) and dynamic reaction mechanisms and are composed of a series of multi-objective optimization procedures located in different environments. The objectives in a MOP are conflicting with each other, so it is essential to explore the trade-off solutions set of all objectives. Different from MOPs, in a DMOP, the objective functions as well as variables maybe change as the time varies from one environment to another. It follows that solving DMOPs is to deal with a series of MOPs generated at different times and obtain the Pareto optimal solutions for each MOP. One optimization procedure is always followed by another and is expected to provide a high-quality Pareto optimal set for each of these static MOPs before the time $t$ changes. As a result, it has become a popular method to divide DMOPs into a sequence of MOPs.

In a DMOP, it is always assumed that the optimization problems at adjacent moments may have certain dependence. Hence, the optimal solutions in the previous moments are always reused to generate the initial population at new environment. In order to make full use of the efficient information of better individuals at past moments, a number of prediction models have been developed. These prediction methods increase the exploitation capability of the population for new solutions in some degree, and by which some classical multi-objective evolutionary algorithms are adopted to deal with DMOPs.

Sahmoud and Topcuoglu ^[10] presented a method for judging change types by calculating the difference of non-dominated solutions in two adjacent environments, and some other similar approaches have also been proposed in the past years ^[11]. These schemes are feasible and efficient when the change of solution sets is regular and stable. However, when the states change sharply, the optimization solutions may deviate significantly from the prediction area. Most of the existing approaches have limitations to deal with this kind DMOPs. This paper intends to develop a strategy combining center and multi-direction that can capture the complex environmental situations.

In our study, we consider the change type of environment with a new method. When the new environment happens, according to the characteristic of the two nearest historical moments, the change type and corresponding different predicting tactics are proposed. The main contributions of this paper are listed as follows:

1) A new method of evaluating individuals is adopted. Rank values of each individual for all objective functions are recorded. Then, solutions are distinguished by their corresponding rank sums, which are used to guide the processes of optimization and dynamic response of the algorithm.

2) On the basis of statistical information of rank sums for all individuals in two adjacent environments, translational change and non-translational change are determined. Different prediction strategies relying on different environment change types are adaptively adopted to respond the new environments.

3) Experiments have been conducted to analyze the performance of our proposed algorithm MOEA/D/CMDS compared with others under different environmental change frequencies and severities.

Next, the related work and preliminaries about DMOPs are introduced. The proposed algorithm is described in Section 4. The test instances, compared algorithms, performance metrics, and parameter settings are presented in Section 5. A series of experiment results and analysis are provided in Section 6. The conclusion and further research are expressed in Section 7.

2.   Related research

In recent years, a large number of researchers proposed dozens of effective methods to solve DMOPs. According to their properties and characteristics, these approaches can be categorized by some representative schemes adopted in algorithm development, such as diversity, memory, prediction, decision variable classification and others. For classical DMOPs, most of these methods are competitive in keeping convergence, or remain superior in decreasing the loss of diversity.

2.1.   Methods based on diversity

In DMOPs, once a new environment occurs, the populations have been convergent. As a result, the optimization algorithm always lacks the exploration ability, and then hardly jumps out of the local region. In order to solve this issue, some techniques, such as random initialization ^[12], hyper mutation ^[13], and immigration strategy ^[14], were used to overcome diversity loss. Deb et al. ^[15] proposed DNSGA-II-A and DNSGA-II-B to increase the population diversity. In DNSGA-II-A, a part of the initial individuals of the population were regenerated to replace the original ones. In DNSGA-II-B, a certain proportion of individuals in the current population was updated by the mutation operation. Woldesenbet and Yen presented a triggered hyper mutation method ^[16]. A new mutation strategy ^[17] was combined with the prediction model to increase the diversity of individuals. In ^[18], Camara and Ortega selected some non-dominated individuals by crowding distance mechanisms, and these individuals were taken as candidates for mutation. Furthermore, a diversity preservation strategy was stated by Ruan et al. ^[19]. Most of these methods are good at keeping diversity, but they seem ineffective for complicated and unpredictable environments.

2.2.   Methods based on convergence

With the purpose of obtaining high quality initial population to speed up convergence, the strategies based on memory and prediction are widely used in DMOPs. Memory-based approaches are available for periodic or cycle environment changes. In this framework, the known historical optimal information of individuals is always utilized to respond to the new environment. Liang et al. ^[20] proposed a hybrid of memory (HMPS) and prediction strategies to improve the performance of the algorithm on the basis of whether the changes were similar to the previous ones. A new prediction and memory strategy (PMS) was applied by Peng et al. ^[21] as an effective response method. In ^[22], Wang and Li considered two memory schemes. The first one selected individuals from the archive and reused them. The second scheme selected individuals from the archive, and then used Gaussian perturbation on these individuals to produce new individuals which survived in the new environment. In ^[23], Koo et al. came up with a novel memory technique by sieving individuals with lesser crowding distance and generating new ones to act as elements for initial population. These schemes, as memory strategies, have been used for solving DMOPs and showed better performance than earlier approaches. Nevertheless, how to take full advantage of existing data is still a key but hard issue. To adopt redundant history information may be computation-expensive and even useless, whereas too little historical information could not provide useful value for exploring new environment.

To fully take advantage of previous experience to guide new search direction, some techniques based on prediction have been utilized generally in DMOPs. A large number of experiment results show that effective prediction mechanism is helpful to accelerate the convergence of the algorithm. Strategies using prediction can learn the change patterns from the previous environment and predict the new locations of the optimal solutions, which can improve the efficiency of the algorithm. In ^[24], PS was conceived to be constructed by a center as well as manifolds (PPS). The known centers were serving to predict new centers, and the subsequent manifolds were estimated by using the previous manifolds. New centers and corresponding manifolds were taken as new individuals for the next generation of the evolution. Zou et al. ^[25] proposed a prediction strategy (CKPS) in which individuals were classified into three categories: non-dominated points, knee points and adaptive random points. The non-dominated ones with existing center point were going to participate in predicting. The evolutionary direction was controlled by knee points. The number of random points was determined by the severity of environmental changes. Kalman filter (KF) strategy was employed to look for PS and PF in the new environment by Muruganantham et al. ^[26]. A multi-directional prediction strategy (MDP) was designed by predicting the chosen representative individuals to guide the evolution direction ^[27]. Cao et al. ^[28] used a differential model to predict the locations of new individuals, in which three adjacent centers were adopted from different environments. In ^[29], Zheng et al. tested the effect of different prediction methods with different decision variables on the convergence and distribution of individuals, respectively. Li et al. ^[30] put forward a dual prediction strategy with an inverse model which predicted individuals in objective space as well as decision space, and good information of objective space was used to provide feedback for decision space. Rambabu et al. ^[31] built a mixture-of-experts-based ensemble framework which combined diversified types of predictors with gating network to describe prediction performance at different environments. In a word, these prediction methods try to get close to the new PS and PF by generating superior individuals to enhance the exploration capability of algorithms. However, the prediction error may be true if the historical information is not enough or predicting model is not inappropriate. Hence, it is necessary to establish a reasonable prediction model to improve the search performance of the algorithm.

2.3.   Methods based on other schemes

Except the methods based on diversity and convergence, there are other strategies to tackle DMOPs. For example, Xie et al. ^[32] adopted decision-variables-classification based cooperative evolution to improve the performance of the algorithm. In DMOPs, the probability distributions of the decision variables may change in different environments. On this basis, Liang et al. ^[33] used the Spearman rank correlation coefficient (SRCC) to measure the correlations between decision variables and objectives functions, which could efficiently produce potential high-quality points. In ^[34], Chen et al. proposed a domain adaptation learning strategy by introducing a mapping matrix between the search spaces of past and current environments. Zhang and Yang ^[35], sampled individuals by classification-based sampling to enhance population diversity and used the method of probability-based space shrinkage to generate high quality solutions. Zhang et al. ^[36] made use of the precision controllable mutation and the simulated isotropic magnetic particles niching to make individuals approximate the entire Pareto front automatically. In addition, other effective strategies have also been proposed in recent years. These methods above can enhance the search performance of algorithms to some extent.

3.   Preliminaries

3.1.   Basic definitions

Definition 1. Time parameter

The time $t$, as a discrete parameter, is defined by the following mathematical formula ^[37]:

where $\tau$ is the number of generations, ${n}_{t}$ means change severity, ${\mathrm{\tau }}_{\mathrm{t}}$ is used to describe the change frequency of environment, and parameters ${n}_{t}$ and ${\tau }_{t}$ control the degree of environmental change.

Definition 2. Dynamic Pareto dominance ^[37]

Suppose $u\left(t\right)$ and $v\left(t\right)$ are arbitrary individuals from decision space $\Omega$ at time $t$. $u\left(t\right)$ is said to dominate $v\left(t\right)$ (denoted by $u\left(t\right)\prec v\left(t\right)$) if and only if the relational expressions ${f}_{i}\left(u, t\right)\le {f}_{i}\left(v, t\right), \forall i\in \left\{1, 2, 3, \cdots, m\right\}, and{f}_{j}\left(u\left(t\right)\right) < {f}_{j}\left(v\left(t\right)\right)$, $\exists j\in \left\{1, 2, 3, \cdots, m\right\}$, are true. $u\left(t\right)$ is called a non-dominated solution if there is not a solution which can dominate u(t).

Definition 3. Pareto optimal solution set

Pareto optimal solutions set at time $t$ is composed of selected non-dominated individuals and denoted by ${PS}_{t}$, which can be expressed by this form:

Definition 4. Pareto optimal front

The Pareto optimal front at time $t$, known as ${PF}_{t}$, can be expressed by the set as below:

Definition 5. The types of DMOPs

According to the dynamic change characteristic of the PS and PF, Farina et al. ^[38] classified DMOPs into four types.

Type I: The PS changes over time but the PF is fixed.

Type II: Both the PS and PF change over time.

Type III: The PS remains fixed, while the PF changes over time.

Type IV: Both the PS and PF remain fixed.

Although the fourth type appears on occasion, the optimal solutions keep constant. As a result, we mainly focus on studying the first three types.

Definition 6. Prediction based on center

Given the stability of the geometrical center, the prediction schemes based on centers are widely used in dynamic optimization algorithms. This is because the moving direction of individuals are consistent with centers in some cases. The center of the PS can be expressed as follows:

where the number of solutions in ${PS}_{t}$ is denoted by $\left|{PS}_{t}\right|$, $x$ is an element pertaining to the set ${PS}_{t}$. This type of prediction model works well for a class of DMOPs in which the change of PS can be described by the center of the PS.

3.2.   Change detection

Dynamic multi-objective optimization always involves change detection, change response and other related procedures. There are two common methods which are used in DMOPs to probe the change of environment. The first one is a fixed detectors approach, which always evaluates a number of detector vectors to observe whether the objective values of them undergo changes. This approach needs to re-evaluate the objective function values at every generation and is not suitable for noisy environments. The second method is a behavior-based approach ^[39,40], which involves the discrepancy and the statistical information of solutions in objective space. The behavior-based detection approach is popular because there are no additional fitness evaluations. As a result, we prefer the behavior-based scheme to detect whether the new environment appears. Once environment changes are detected, most of the strategies described above based on diversities and convergence belong to common change response methods.

4.   Proposed algorithm

In the section, a new method of dividing change type is proposed. The change types of environments are described in Figure 1. They are divided into translational changes and non-translational changes. Predicting models mainly depend on different change types. Therefore, for the translational change, a predicting strategy with center is adopted. Conversely, the multi-direction strategy is used to generate new individuals adaptively. The predicting model is named CMDS. In our studies, MOEA/D is used to execute evolution and the Tchebycheff approach acts as a decomposition method in static optimization stage ^[41]. Figure 2 shows the framework of the proposed algorithm (MOEA/D/CMDS). Next, we describe the detailed process of the algorithm. The whole structure of MOEA/D/CMDS is shown in Algorithm 1.

4.1.   Types of changes

For each individual ${x}_{i}^{t}, i\in \left\{1, 2, \cdots, N\right\}$ at time $t$, its corresponding objective function value is ${f}_{k}\left({x}_{i}^{t}\right), k\in \left\{\mathrm{1, 2}, \cdots, m\right\}$. All elements in set $\left\{{f}_{k}\left({x}_{1}^{t}\right), {f}_{k}\left({x}_{2}^{t}\right), \cdots, {f}_{k}\left({x}_{N}^{t}\right)\right\}$ are sorted from the smallest value to the largest one. We record the sequential number of ${x}_{i}^{t}$ with ${R}_{i}^{kt}$. In other words, the rank value of ${x}_{i}^{t}$ is ${R}_{i}^{kt}$ for the k-th objective at moment $t$. The smaller sequential number means smaller function value for each solution on the corresponding objective. The sum of ${R}_{i}^{kt}$ for individual $i$ in each objective function is denoted by ${R}_{i}^{t}$, and its expression is as follows:

where the maximum value of ${R}_{i}^{t}$ is ${max}^{t}$ and the minimum value of it is ${min}^{t}$ for the varying $i$. The sequence ${{R}^{t} = (R}_{1}^{t}, {R}_{2}^{t}, \cdots, {R}_{N}^{t})$ represents the rank values of all solutions at moment $t$. As it should be, the sequence ${{R}^{t-1} = (R}_{1}^{t-1}, {R}_{2}^{t-1}, \cdots, {R}_{N}^{t-1})$ displays the rank values at moment $t-1$. The intensity variation of the rank values corresponding to individual $i$ is represented by the difference between ${R}_{i}^{t}$ and ${R}_{i}^{t-1}$. That is to say, the location variation degree for the i-th individual at $t$ moment and $t-1$ moment is measured by the value $\left|{R}_{i}^{t}-{R}_{i}^{t-1}\right|$. When the absolute value is small enough, we can assume that the change of location is not obvious for corresponding solution, or else the change of location is sharp. In order to distinguish the change type, the average value of $\left|{R}_{i}^{t}-{R}_{i}^{t-1}\right|(i = 1, 2, \cdots, N)$ is calculated and denoted by ${V(R}_{t})$. It is expressed below:

we construct threshold value ${\mathrm{\alpha }}_{\mathrm{t}}$ which is below at time $t$ :

where ${\alpha }_{t}$ varies with environment. At the same time, as the dimensionality $m$ of the objective function increases, ${max}^{t}$ and ${min}^{t}$ grow intensely. Hence, the difference of rank sums between the best solutions and the worst solutions is used to depict a change degree for extreme individuals. Smaller value of it means better stability of population relative to the increment of objective function number $m$. If the position movement of all individuals between moments $t-1$ and $t$ is not drastic, ${V(R}_{t})$ must be smaller. Conversely, the larger ${V(R}_{t})$ is true. According to the above idea, the threshold value ${\alpha }_{t}$ can be set. When ${V(R}_{t})$ is less than ${\alpha }_{t}$, the change type of environment is translational; Otherwise, it is considered as non-translational. For example, in the Figure 1, the black dots represent individuals at $t-1$ moment, and the red dots show solutions at $t$ moment. The solutions are noted as ${x}_{1, }{x}_{2, }{x}_{3, }{x}_{4, }{x}_{5}$ from left to right at moment. Arrows indicate the direction of individual evolution. In part (a) of Figure 1, rank sums of each solution do not vary. Hence, ${V(R}_{t})$ is smaller than ${\alpha }_{t}$ and the change type of environment is translational. However, in part (b) of Figure 1, individual ${x}_{2}$ has the maximum rank sum $8$ and individual ${x}_{1}$ has the minimum rank sum $3$ at the current moment $t$. It is obvious that ${V(R}_{t})$ is larger than ${\alpha }_{t}$, so the change type of environment is determined. The frame of identifying the change type of new environment is stated in Algorithm 2.

4.2.   Prediction strategies

The intention of prediction is to make the individuals quickly converge to the new PF before the next environment change is detected. Effective prediction strategies can reduce computation costs and promote algorithm convergence. Cao et al. proposed a second-order difference model to predict the moving trends of centroid of the solutions obtained ^[28], but the distinctions of individuals were neglected. The centers adopted by them were determined with all ones of the population. In our paper, solutions are selected, making up for the above shortcomings, and we employ a new prediction strategy which responds to the environment in terms of environment change type.

4.3.   Response to translational change

Due to the high efficiency of using center predicting, when the change is translational, there are no extreme individuals whose rank sums alter intensely. Thus, we can conclude that the present environment could have great similarity with previous one. The change of environment is regular and not sharp. That is to say, the change is relatively stable. Hence, we can assume that the movement direction of the individual is consistent with that of the center of PS.

In this paper, a new estimation method of PS is put forward. For the sequence ${R}^{t}$, if the value of ${R}_{i}^{t}$ is relatively small, it follows that individual $i$ has good characteristic. All individuals with the same properties as individual $i$ constitute the set ${B}_{t}$, which approximates ${PS}_{t}$ at the moment $t$. The individuals which do not belong to set ${B}_{t}$ are deemed to be in the set ${A}_{t}$. It is apparent that present population is constituted by the elements of the set ${A}_{t}$ and ${B}_{t}$. How to determine elements in the set ${B}_{t}$? The smaller rank value of individual implies the better proximity to PS. However, if there are too few elements in set ${B}_{t}$, the characteristic of the population could not be reflected fully. For the same reason, too many elements do not exactly reveal the superiority of individuals with lower rank values for ${B}_{t}$. Taking these factors into consideration, the threshold value ${\beta }_{t}$ is calculated in the following way:

It is evident that ${\beta }_{t}$ varies over time $t$, and the average allocation of the difference between ${max}^{t}$ and ${min}^{t}$ for each objective function is used to balance discrepancies of whole objectives. If the value of ${R}_{i}^{t}$ less than ${\beta }_{t}$ is true, corresponding individual $i$ is put into set ${B}_{t}$. By the same principle, the sets ${B}_{t-1}$ and ${B}_{t-2}$ can be determined by the corresponding ${\beta }_{t-1}$ and ${\beta }_{t-2}$ which are threshold values corresponding to moments $t-1$ and $t-2$. Apparently, individuals with higher quality are gathered in sets ${B}_{t}$, ${B}_{t-1}$ and ${B}_{t-2}$. Then, the centers of these sets are separately denoted by ${C(B}_{t})$, ${C(B}_{t-1}),$ and ${C(B}_{t-2})$. Without loss of generality, in order to save computation cost, the previous three consecutive centers are used to construct new movement direction. In the process of predicting, maybe prediction error exists. To overcome this shortcoming, Gaussian perturbation with mean value $0$ and standard deviation $\delta \left(t\right)$ is introduced. The prediction formula is designed as follows:

where ${r}_{1}$ is a parameter randomly taken in $\left(\mathrm{0, 1}\right)$ and $\delta \left(t\right)$ is expressed in the following form:

It is obvious that smaller value ${max}^{t}-{min}^{t}$ means more clustered individuals, so the disturbance amplitude should be larger. This symbolizes that larger value of $\delta \left(t\right)$ is established. Otherwise, the distribution of individuals is relatively scattered. The population has relatively better diversities, then only smaller disturbance amplitude is required. In a word, Gaussian perturbation is introduced to improve the search and exploration capability of population. In this response mechanism, the present population is divided into two subpopulations. Solutions of the first subpopulation are used to execute prediction. The individuals in the second subpopulation serve as the initial ones in the new environment directly. It is worth noting that there is no stored historical information to meet with the prediction formula of Eq (10) when the algorithm runs at the earliest three moments. In this case, all of the individuals are randomly initialized to respond the new environment change.

4.4.   Response to non-translational change

When the non-translational change is detected, the population is divided into three parts. For different parts, the corresponding individual generating scheme is adopted. The set ${B}_{t}$ is as described above at present moment. It is easy to find that subpopulations ${B}_{t}$ and ${A}_{t}$ are distinguished by the value of ${\beta }_{t}$. For a given individual ${x}_{u}^{t}$ in set ${B}_{t}$, considering each individual ${x}_{v}^{t-1}$, $v = 1, 2, 3, \cdots, \left|{B}_{t-1}\right|$ in set ${B}_{t-1}$, the Euclidean distance ${dist(x}_{u}^{t}, {x}_{v}^{t-1})$ between ${x}_{u}^{t}$ and ${x}_{v}^{t-1}$ is calculated. The individual ${x}_{v}^{t-1}$ with minimum ${dist(x}_{u}^{t}, {x}_{v}^{t-1})$ is recorded as ${x}_{{u}^{*}}^{t-1}$. It follows that ${\mathrm{x}}_{{u}^{*}}^{\mathrm{t}-1}$ is the nearest point to ${x}_{u}^{t}$ for all individuals in set ${B}_{t-1}$. Using the same way, the individual ${x}_{{u}^{*}}^{t-2}$ could be found in the set ${B}_{t-2}$. The Euclidean distance of ${x}_{u}^{t}$ and ${x}_{{u}^{*}}^{t-2}$ is the minimum for all solutions in the set ${B}_{t-2}$. Next, a new prediction model which is called multi-direction strategy is proposed for each individual ${x}_{u}^{t}$ in set ${B}_{t}$.

where ${r}_{2}$ is a random number from $\left(\mathrm{0, 1}\right)$. We utilize pointwise prediction model for each element in the set ${B}_{t}$. Each individual in ${B}_{t}$ has its direction. Obviously, the individuals which do not take part in this prediction belong to the set ${A}_{t}$. We consider the individuals whose rank sum is ${min}^{t}$ at present time $t$. In many cases, more than one individual has such feature, so these ones are known as the best solutions, which constitute the set ${G}_{t}$. We assume that the center of ${G}_{t}$ is $C\left({G}_{t}\right)$. Then, we categorize the elements in the set ${A}_{t}$ according to random probability to increase the diversity of the population. For individuals in set ${A}_{t}$, the first part solutions which are regarded as belonging to set ${A}_{t1}$ are selected by a random probability to accept the below prediction strategy:

where ${C}_{tt}$ is the center of the whole population and $\gamma$ is a random number from $\left(\mathrm{0, 1}\right)$. The second part solutions which constitute the set ${A}_{t2}$ are the rest ones in the set ${A}_{t}$. The way of generating new individuals for it is random initialization. The method is showed as below:

where ${Lb}_{i}$ and ${Ub}_{i}$ which are given in the experiment are the upper bound and lower bound of each dimension variable, $i = 1, 2, \cdots, n.$ $i$ represents dimension of individuals ${x}_{ri}^{t}$ and ${x}_{ri}^{t+1}$. $\theta$ is a random number relying on $[0, \mathrm{ }1]$. The detailed description about response strategies for two change types of environments of CMDS is described in Algorithm 3.

In general, a method relying on rank sums of individual is designed in this paper. Different prediction strategies are designed for individuals of different categories. Next, the results of the experiment are showed in the following section.

5.   Experiment design

5.1.   Test instances

The proposed algorithm MOEAD/D/CMDS is tested on 20 benchmark instances, including FDA1–FDA5 test suites ^[38], dMOP1–dMOP3 ^[42] test suites, ZJZ test suites (F5–F8) and JY1–JY8 test suites ^[24]. By extending FDA problems, the dMOP test suites are acquired. A linear relationship exists among the decision variables of FDA and dMOP. However, there are nonlinear correlations in decision variables of F5–F8. For JY test suites, there are mixed PF and intricate relationships among decision variables which are all challenging features for dynamic multi-objectives optimization.

5.2.   Compared algorithms

The proposed algorithm is compared with MOEA/D/DM ^[28], MOEA/D ^[41], MOEA/D/KF ^[26], DNSGA-II-A ^[15] and DSS ^[43]. MOEA/D/DM used a second-order center difference model to predict the moving trends of part individuals. Classical multi-objective optimization algorithm MOEA/D was employed for DMOPs after being slightly modified, wherein a part of individuals was randomly initialized. As is known to all, the same idea was applied to DNSGA-II-A. In the algorithm MOEA/D/KF, Kalman filter model was used to search the new direction of solutions. However, different from the above DMOPs algorithms, a directed local search method with linear prediction model is taken into DSS.

5.3.   Performance metrics

Performance metrics can evaluate convergence, distribution, and diversity of the obtained population. The generation distance $\left(GD\right)$ has been widely used in evaluating algorithm and was introduced in ^[42]. The expression form of $GD$ is given as follows:

where the ${PF}_{t}$ is consisted by evenly distributed points of the true PF at $t$ moment and ${P}_{t}$ is a set which is approximate to ${PF}_{t}$. The definition of $d\left({PF}_{t, }v\right)$ is as follows:

where $d\left({PF}_{t, }v\right)$ is the smallest Euclidean distance between individual $v$ and the set ${PF}_{t}$. $\left|\bullet \right|$ represents cardinality of the corresponding set. The smaller $GD$ of the algorithm symbolizes the better convergence. The $IGD$ ^[19,44] is a modified version of $GD$. The expression of it is as below:

and $d\left(u, {P}_{t}\right)$ can be expressed as:

$d\left(u, {P}_{t}\right)$ means the minimum Euclidean distance between point $u$ on the ${PF}_{t}$ and the set ${P}_{t}$. As is known to all, $IGD$ is a useful index to measure the convergence and diversity of the obtained solutions.

The average value of $IGD$ is defined as $MIGD$ metric for all time windows in a single run ^[45].

where $T$ is a set and represents all moments. Its cardinality is $\left|T\right|$. The smaller value of $MIGD$ means the better performance of the algorithm. In this paper, we employ $MIGD$ value to measure the performance of the algorithm MOEA/D/CMDS.

5.4.   Parameter setting

Considering the change frequency and change severity affect the performance of algorithm. The experiment is tested on ${n}_{t} = 5, 10, {\tau }_{t} = 5, 10.$ Hence, there are $4$ types for each problem. In order to maintain fairness, the parameters in the comparison algorithms are taken from original references. The size of the population is $100$ and 300 for bi-objective problems and tri-objective problems. The values of parameters ${r}_{1}$ and ${r}_{2}$ are equal to 0.6. At the same time, the value of $\gamma$ is set to 0.8. Relatively larger values of ${r}_{1}$, ${r}_{2}$, and $\gamma$ reflect that present environment data is dominant for generating new solutions. ${Lb}_{1}$ and ${Ub}_{1}$ are taken as 0 and 1, separately. Meanwhile, ${Lb}_{i}$ is assigned as $-1$ and ${Ub}_{i}$ is set as 1 when $i$ changes from 2 to $n$.

In DE operator, $CR$ is fixed as 0.5 and $F$ is set to 0.5. In the polynomial mutation operation, the value of parameter $\eta$ is 20. ${p}_{m} = 1∕n$ is true, in which $n$ represents the dimensions of the decision space. The size of the neighborhood is 20. The coefficient $\delta$ is 0.8. The size of ${n}_{r}$ is the same as the cardinality of mating pool E, where parent individuals are selected to generate new ones. The total number of iterations is set to ${40\tau }_{t}$, which assure 40 environmental changes to be carried out in each run. The proposed algorithm is executed independently 30 times for each example, and the $MIGD$ value is taken as a measure metric, the smaller $MIGD$, the better performance of MOEA/D/CMDS.

6.   Experiment result analysis

The Tables 1–4 show the experiment results of each type of problem on $MIGD$ values and their standard deviations. The best results are marked in bold. The detailed analysis about the test result is as follows. The curves of tracking the $IGD$ trend over 40 environment changes for 20 runs are showed in the Figures 3–5.

6.1.   Results on FDA and dMOP problems

This section shows the experiment results about test functions FDA and dMOP. The PS of FDA1 changes as a sine curve and the PF keeps fixed. MOEA/D/CMDS performs better than or close to MOEA/D/DE except the case ${n}_{t} = 5, {\tau }_{t} = 10$. For FDA2 problem, the shape of PF is convex and non-convex alternation in the objective space. At this case, MOEA/D/CMDS is slightly superior to or equivalent to other compared algorithms for all cases, which indicates that prediction strategy can improve the search ability of the algorithm. In FDA3 and FDA5, different density distributions exist in distinct segments of the PF. These two types of changes don't follow any assumed models, which cause certain challenges for all algorithms. However, the individuals with smaller rank sums play a positive guidance in the former adjacent environment for high severity change, which show the effectiveness of MOEA/D/CMDS. The test function FDA4 is a tri-objective problem whose PF and PS change with time. In this case, the moving direction provided by MOEA/D/CMDS is a disadvantage because the change of position sometimes inaccurately describes the movement direction of individuals. The results of FDA4 for all parameters are inferior to MOEA/D/DM, but they are obviously better than other compared algorithms. Specific experimental results are shown by Table 1.

The statistical results about dMOP are listed in Table 2. The change type of dMOP1 belongs to Type III which keeps static PS, whereas the PF and PS of dMOP2 are all time-varying. For these two instances, when the change frequency parameter ${\tau }_{t}$ is equal to 5, MOEA/D/CMDS performs well because of powerful guidance of better history solutions. At the same time, if the value of ${\tau }_{t}$ is equal to 10, the degree of change of individual position is very small. Hence, the loss of diversity is obvious, which causes the degradation of algorithm performance. In dMOP3, only PS changes over time, and the variables which control the spread of the PF vary randomly. This mechanism of change brings certain difficulties for prediction schemes. Good performance of MOEA/D/CMDS shows its superiority in dMOP3. Figure 3 presents the trajectories of IGD on these two test problems.

6.2.   Results on ZJZ problems

The ZJZ test instances include F5–F8. The nonlinear relation is obvious between decision variables. Besides, according to the value of ${V(R}_{t})$, the changes type of PF is non-translational. The change of the position for individuals is sharp in two adjacent environments. Hence, it is difficult to fully utilize the information of the previous time windows. Therefore, ZJZ problems cause more challenges for DMOPs, especially when they face higher frequency and more severity changes. From Table 3, the algorithm proposed in this paper is very advantageous compared with other rivals except F8 which is a test function of three objectives. It is obvious that the $\mathrm{M}\mathrm{I}\mathrm{G}\mathrm{D}$ values of MOEA/D/CMDS are far less than competitors. But for F8, the individuals which have smaller rank sums in the previous environment can't guide new search direction well, so MOEA/D/CMDS is disadvantage. In a word, MOEA/D/CMDS performs well on most of the test instances for ZJZ problems, and the value of $\mathrm{M}\mathrm{I}\mathrm{G}\mathrm{D}$ has superiority. The $\mathrm{I}\mathrm{G}\mathrm{D}$ curves of tracking characteristic are showed in Figure 4, the results of experiment are provided by Table 3.

6.3.   Results on JY problems

The results of the JY test are showed in Table 4. MOEA/D/CMDS has the better superiority on four problems JY3, JY4, JY6, and JY7, but it performs a little worse than the rest test problems. The decision variables of JY1 and JY2 are irrelevant and the center difference model of the whole population in MOEA/D/DM could improve the performance of the algorithm. For the JY3 problem, any two decision variables have time-varying non-monotonic dependencies. As time goes on, the PS becomes more and more complicated. At this case, MOEA/D/CMDS considers the nearest ones in B_(t-1) and B_(t-2) for each individual in B_t in order to make full use of history information. Of course, individuals with poor quality are discarded at the same time. Therefore, solutions from previous moments with smaller rank sums could improve the exploration ability of algorithm. The MIGD values of JY3 are far less than other compared algorithms for all cases. The PF of JY4 is discontinuous and its disconnected parts change with time. It is difficult to cover all the PF components. Although this type of problem is difficult, the method adopted in this paper is competitive when the environment changes severely. In other cases, the results are similar to MOEA/D/DM because the PF of JY5 keeps stationary. Thus, MOEA/D/CMDS performs poorly on this kind of problem. For the test functions JY6 and JY7, they are multimodal problems and the PF shapes of them appear concave and convex alternately. The number of local optima changes for JY6, but the opposite situation is true for JY7. Local optimum increases more evaluations to acquire the PS when the environment does not change. As it should be, algorithms may fall into local optimum, which affects their performance. In short, dealing with these two types of problems is very difficult. The PS of JY8 is static, and yet the segments of PF change with time. When the change of environment is fast and severe, the experiment results are approximate to the results of MOEA/D/DM. In other cases, MOEA/D/CMDS performs worse than MOEA/D/DM and MOEA/D, but it is obviously superior to most of the comparison algorithms. Features of JY instances are reflected by Figure 5. The experiment results about JY series are shown by Table 4.

We mark the computational results with "+", "-" or "≈", which indicate MOEA/D/CMDS is better than, worse than or equal to the compared algorithms. In a word, on FDA test set, the proposed algorithm shows its superiority on 10 among 15 test problems. For dMOP and ZJZ test problems, the proposed algorithm defeats other compared algorithms on at least three fourths of the problems. For the test functions JY, half test results of MOEA/D/CMDS are evidently superior to other compared approaches.

7.   Conclusions and future work

In this paper, a new prediction scheme named CMDS is first presented to adaptively select prediction strategy to respond new environment. When environment change happens, CMDS method firstly judges the type of change, which is translational or no-translational. Next, the way of generating individual is determined in order to provide initial individuals at next moment. As a result, CMDS method is integrated into the frame of MOEA/D. From the experiment results, MOEA/D/CMDS can show its superiority in tracking the changing PS or PF.

In the future research, we will be devoted to studying more efficient prediction mechanisms to handle more complex dynamic optimization problems, and further extend to large-scale practical problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant (No.61966030) and the Qinghai Provincial Key Laboratory of Internet of Things (No.2022-ZJ-Y21).

Conflict of interest

The authors declare that there are no conflicting interests.