Robust Optimization Over Time by Estimating Robustness of Promising Regions

Many real-world optimization problems are dynamic. The field of robust optimization over time (ROOT) deals with dynamic optimization problems in which frequent changes of the deployed solution are undesirable. This can be due to the high cost of switching the deployed solutions, the limitation of the needed resources to deploy such new solutions, and/or the system being intolerant toward frequent changes of the deployed solution. In the considered ROOT problems in this article, the main goal is to find solutions that maximize the average number of environments where they remain acceptable. In the state-of-the-art methods developed to tackle these problems, the decision makers/metrics used to select solutions for deployment mostly make simplifying assumptions about the problem instances. Besides, the current methods all use the population control components, which have been originally designed for tracking the global optimum over time without taking any robustness considerations into account. In this article, a multipopulation ROOT method is proposed with two novel components: 1) a robustness estimation component that estimates robustness of the promising regions and 2) a dual-mode computational resource allocation component to manage subpopulations by taking several factors, including robustness, into account. Our experimental results demonstrate the superiority of the proposed method over other state-of-the-art approaches.


I. INTRODUCTION
S EARCH spaces of many real-world optimization problems are dynamic. To tackle optimization problems in dynamic environments, it is important that the optimization algorithms can efficiently react to the environment changes [1], [2]. A dynamic optimization problem (DOP) can be defined as where f is the objective function, t ∈ [0, T] is the time index, x is a solution in the search space, and α is a vector of time-dependent control parameters of the objective function. Almost all existing works in the field of DOPs consider problems whose environmental changes occur only at discrete-time steps, i.e., t ∈ {1, . . . , T}. For a DOP with T environmental states, there is a sequence of T stationary environments f x, α (k) T k=1 = f x, α (1) , f x, α (2) , . . . , f x, α (T) . ( Some main characteristics of a DOP are defined based on its change severity and frequency. Hence, these two are among the main criteria used for classifying DOPs [2]. The change severity shows how much the morphology of the problem space changes after each environmental change. In the DOP literature, it is mostly assumed that environmental changes are not highly severe and there is a degree of similarity between consecutive environmental states. This is the case for many practical applications [3]- [5]. In real-world DOPs, change frequency is defined based on the duration of the time interval between environmental changes which depends on the nature of the events that cause the environmental changes. In some problems, the duration of this time interval can be very short, which results in higher change frequency (e.g., the short time gaps between changing demands/customers in some dynamic covering location problems [6]). Change frequency is usually very low (i.e., longer time gap between changes) for problems where environmental changes are caused by accidents or faults in parts of the system. Evolutionary algorithms are commonly used for tackling DOPs [1], [5], [7]. However, since these algorithms are originally designed for solving static optimization problems, they cannot directly be used for DOPs. This is due to the challenges caused by the environmental changes in DOPs, which are: global and local diversity loss, outdated stored fitness values, 1 and limited number of fitness evaluations that can be performed in each environment (i.e., between two consecutive environmental changes) 2 [8]. In order to avoid influences of some factors, such as hardware, compiler, and programming skills, the academic community of evolutionary dynamic optimization (EDO) uses the number of fitness evaluations as the unit of change frequency instead of time in the real-world applications. In many real-world DOPs, the fitness evaluation is time consuming (e.g., in large-scale problems [9], or simulation optimization). Consequently, a limited number of fitness evaluations can usually be performed in each environment. To address the aforementioned challenges of DOPs, EDO algorithms are often created by augmenting evolutionary algorithms with other algorithmic components to address the DOPs' challenges stated above [2].
The majority of the existing EDOs tackle DOPs by tracking the moving global optimum after each environmental change [3], [10]. However, tracking is impractical for solving many real-world DOPs because frequent change of the deployed solution is not possible. This can be due to different reasons, such as high switching cost [11], [12], limited available resources to deploy new solutions, or system intolerance for frequent changes in the deployed solution [13], [14].
To solve this type of DOPs, Yu et al. [13] proposed an approach called robust optimization over time (ROOT). In ROOT, in order to reduce the number of times that the deployed solution is changed, it is reused (i.e., kept deployed) until its quality degrades to an unacceptable level. Therefore, although a deployed solution in an environment is not necessarily the best solution, it must satisfy a quality-based constraint [15]. In this type of ROOT problems, a new solution must be chosen for deployment when the current deployed solution is no longer acceptable after the last environmental change [10]. We wish to maximize the average number of environments where the previously deployed solutions can be reused and kept deployed [15], [16]. Later, two other types of ROOT problems were investigated: 1) time window-based [17] and 2) switching cost-based [12] ROOT problems. In this article, we focus on the first type, i.e., ROOT problems with a quality-based constraint [13]. Unless otherwise stated, we use the term ROOT to refer to this specific type of ROOT throughout this article.
To solve a ROOT problem, not only should an EDO be capable of addressing the challenges of reacting/responding to the environmental changes but also the challenges of estimating 1 Also called the outdated memory issue in the DOP literature. 2 Also called limited available computational resources in each environment. the robustness and acceptability of solutions in the forthcoming environments. A desirable solution in ROOT can remain acceptable for a higher number of environments in the future. Nevertheless, accurate estimation/prediction of the future acceptability of solutions is very challenging. Depending on how the existing ROOT methods deal with this challenge, they can be categorized into fitness prediction [15], [17] and promising regions' reliability-based approaches [10], [18]. In prediction-based methods, the actual fitness function is altered with a substitute fitness function that considers the predicted fitness values of the candidate solutions in the upcoming environments [15], [17]. It is shown in [18] that using the predicted future fitness values to find robust solutions is too error prone for problem instances generated by the moving peaks benchmark (MPB) [19]. Reliability-based methods, on the other hand, choose the solutions for deployment based on the estimated behavior of the promising regions instead of predicting the future fitness values of solutions [18], [20]. In these methods, some reliable promising regions are determined based on the gathered information by a multipopulation method. Thereafter, a solution is chosen for deployment from a reliable promising region based on a strategy, such as picking the best found solution in the reliable promising region with highest fitness 3 [20], or the best found solution in the promising region with the smallest estimated shift severity [18].
Despite the importance of ROOT in many real-world applications and more than a decade from the first time it was introduced [13], very little attention has been given to the field. A major weakness of the existing methods is that they are all tailored for very simple problems. For example, some of their components are specifically developed for low dimensions [16], simple dynamics [18], and/or regular/smooth search space/promising regions [15], [16], [18]. This creates a gap between the academic research and the real-world applications. Moreover, all existing ROOT methods use EDOs [10], [15], [16], which are originally designed for performing tracking the moving global optimum where robustness is not considered in their population management and control components. Despite the significant role of EDOs in the ROOT methods, little attention has been given to design EDO components which are ROOT-specific.
In this article, we propose two components that address the aforementioned shortcomings: 1) a robustness estimation component and 2) a systematic dual-mode computational resource allocation (CRA) component. By combining the proposed components with those of a multipopulation EDO capable of tracking multiple moving promising regions, a new ROOT method is formed. Thanks to the proposed robustness estimation component, the new ROOT method is not dependent on some oversimplified characteristics of the benchmark problems. This component keeps and transfers the historical knowledge about the covered promising regions and estimates their robustness degrees accordingly. The estimated robustness degree of each promising region is calculated based on the number of environments that previous promising region's summit positions could be reused/kept as the deployed solution until the current environment. The estimated robustness values of the promising regions are used for choosing the next solution for deployment and controlling the subpopulations. A systematic dual-mode robustness-based CRA component for managing the subpopulations is also proposed. Using this component, the proposed method controls the subpopulations based on the estimated robustness of the covered promising regions and the system status. 4 The organization of the remainder of this article is as follows. Section II covers the related works. The proposed method is described in Section III. Section IV explains the experiment setup, including the used benchmark, performance indicator, comparison algorithms, and parameter settings, and also reports experimental results, comparisons, and analysis. Finally, Section V concludes this article.

II. RELATED WORK
A ROOT method is usually constructed by assembling an EDO and some additional components, such as a decision maker [20] or a transformation of the objective function [15], [17]. The literature of EDOs is not covered in this article due to space limitations and the readers can refer to the two-part survey in [2] and [5]. Instead, we only focus on reviewing ROOT methods and the components specifically designed for them.

A. ROOT Problems
Three types of ROOT problems have been investigated in the literature: 1) ROOT problems based on acceptability of the deployed solution in which the deployed solution is kept as long as it remains acceptable [13]; 2) ROOT problems with time window in which the deployed solution is kept during each time window [17]; and 3) ROOT problems based on acceptability of the deployed solution and switching cost in which the deployed solution is kept until it becomes unacceptable or another solution is found whose fitness is considerably higher than that of the deployed solution, making the benefit of switching outweigh the cost [12]. As stated before, in this article, we focus on the first type of ROOT problems where the main objective is to minimize the number of times when the deployed solution is changed.
Given a DOP f (t) ( x) with T environments, the aim of ROOT is to find a set of deployed solutions S = { s 1 , s 2 , . . . , s l } where 1 ≤ l ≤ T. The main objective in ROOT is to minimize l. A deployed solution s i ∈ S is considered robust if it remains acceptable across more than one environment. For a deployed solution s i to be acceptable in the tth environment, f (t) ( s i ) must be greater than a predefined threshold μ [17]. Otherwise, a new solution must be chosen for deployment. This acceptability evaluation approach has been used to determine the acceptability of solutions in the majority of the works in the literature [11], [18], [20], [21].

B. ROOT Methods
Despite the importance of ROOT in tackling many realworld DOPs, this field has not received much attention so far. To the best of our knowledge, there are only three main ROOT methods, which are proposed by Jin et al. [15], Fu et al. [17], and Yazdani et al. [18]. The rest of the works in the field are designed based on these three works [12], [22], [23].
Jin et al. [15] proposed the first ROOT method in 2013. This method uses the predicted fitness values of solutions in a predefined number of future environments to estimate their robustness. The main components of this method include a single-population EDO, a database, an approximator, and a predictor. The single population EDO is responsible for gathering data to train the predictor and also the optimization process. The historical data are archived in the database over time, which is used to train the approximation and prediction methods. Jin's method uses a substitute objective function, which is the accumulation of the actual fitness, predicted, and approximated values.
Fu et al. [17] proposed a survival time-based ROOT method whose main components are the same as Jin's method, however, the used substitute objective function is different where f (t+i) is the prediction function that predicts the future fitness value of x in the (t+i)th environment. According to (3), if the fitness value of a solution in the current environment is less than μ, this solution is considered as a nonrobust solution. Otherwise, robustness value of this solution will be the number of successive future environments in which its fitness values are predicted to remain above μ.
Several multiobjective ROOT methods have been proposed to find robust solutions based on the substitute objective function in (3) and an additional objective function such as average fitness of the deployed solutions or switching cost. In [24], the algorithm tries to find a Pareto front based on the substitute objective functions of average fitness over a predefined time window and survival time. In [11], a multiobjective approach, called ROOT/SC, is proposed for optimizing two objectives, including maximizing the survival time (3) and minimizing the switching cost. The switching cost is defined as the Euclidean distance between the current deployed solution and a candidate solution. ROOT/SC was modified in [21] by adding the current fitness of candidate solutions as the third objective. Another group of multiobjective ROOT methods is designed to find robust Pareto optimal solutions [25]- [27]. In these methods, the main goal is to find Pareto optimal solutions whose performance is acceptable for the current and upcoming environments.
Yazdani et al. [18], [20] proposed a reliability-based ROOT method. Unlike the previous ROOT methods that search for robust solutions based on a substitute objective function, this method works on the search space constructed by the original objective function. The components include a multipopulation EDO capable of locating and tracking multiple moving promising regions and a decision maker that chooses the solutions for deployment. The multipopulation EDO is responsible to gather information about the promising regions (e.g., shift severity and fitness fluctuation degree). Based on the gathered information, the reliability of each covered promising region is determined for choosing the next solution for deployment. To this end, first, the fitness fluctuation degree of each region covered by the ith subpopulation (pop i ) is calculated after each environmental change as follows: is the fitness fluctuation in the tth environmental change, and g * (t−1) i is the best found position by pop i in the t − 1th environment. Thereafter, the average values of τ i in the past environments (τ i ) are used for determining the reliability of the region using If the best found position g * (t) i is reliable, it means that it is expected that its worst possible future fitness value will remain above the threshold μ for at least another environment. In the tth environment, if the previous deployed solution becomes unacceptable, the best found positions in the reliable promising regions (i.e., ρ i = 1) are considered as a set of candidate solutions C for the next deployed solution. In [20], the solution j in C is picked for deployment using where the solution from C that has the highest worst estimated future fitness value is chosen for deployment. Another strategy, used in [18], to choose a solution for deployment is as follows: where s i and h i are the estimated shift and height severity values of the promising region covered by the ith subpopulation, respectively, and s max and h max are the largest s and h values among reliable promising regions, respectively. A major shortcoming of the reviewed methods in this section is that they are all tailored for very simple problems. For example, some of their components are designed for low numbers of dimensions, regular/smooth search space/promising regions, and/or simple dynamics. As described before, both Jin's and Fu's methods use approximation and prediction components to estimate the solutions' future fitness values. On the one hand, it is shown in [10] that using such approximation and prediction methods to estimate the future fitness values of solutions can be error prone. On the other hand, Yazdani's method depends on the accuracy of the estimated fitness fluctuations in (4). Until now, this method has only been tested on MPB whose promising regions are regular/smooth, without ill-conditioning, fully separable, and symmetric peaks [28], with fixed peak relocation length over time and homogeneous dynamics [2]. However, our investigations indicate that the effectiveness of determining reliability of the promising if a solution need to be deployed then 11 Choose the best found solution in the promising regions with the highest estimated robustness;

until stopping criterion is met;
regions based on the estimated fitness fluctuations deteriorates where the problems are more challenging. The last issue concerns the EDOs used in every ROOT methods to perform the optimization in dynamic environments. All existing ROOT methods use EDOs, which are originally designed for tracking the moving global optimum (Jin's and Fu's methods use a single-population EDO presented in [29], and Yazdani's method uses the multipopulation EDO from [30]). Despite the significant role of EDOs in the ROOT methods, little attention has been given into designing some components of EDOs, which take ROOT's considerations into account.

III. ROBUSTNESS ESTIMATION AND COMPUTATIONAL RESOURCE ALLOCATION FOR ROOT
In this section, to address the shortcomings of the existing ROOT methods, we propose two new components.
1) A robustness estimation component that uses an explicit archive to keep and transfer historical knowledge about the covered promising regions for estimating their robustness degrees. The estimated robustness values of the promising regions are used for choosing the next solution for deployment and controlling the subpopulations. 2) A systematic dual-mode robustness-based CRA component that picks the subpopulations to run in each iteration in order to manage the consumption of the fitness evaluations. This is done according to several factors, including the system status, the estimated robustness of each covered promising region, and roles and task achievements of subpopulations. We integrate the proposed components into a state-of-theart multipopulation EDO to construct a ROOT algorithm. Algorithm 1 shows a high level procedure of the resulting ROOT method. As shown in this pseudocode, the proposed robustness estimation component is executed after each environmental change (line 8). The proposed resource allocation is executed at the beginning of each iteration (line 3) to pick the subpopulations to be run in the current iteration.
The high-level components of a multipopulation EDO can be classified into change reaction (line 9), optimization (line 5), and population management and diversity controlling components (line 6), which usually include the mechanisms used for dividing the population into subpopulations, removing/randomizing redundant individuals/subpopulations, generating new subpopulations, and increasing/maintaining global diversity. 5 Finally, a decision maker is used (line 11) to choose a solution for deployment.
In the rest of this section, we describe the following. 1) Suitable multipopulation EDOs that can be equipped with the proposed components to form ROOT algorithms.
2) The proposed robustness estimation component.

A. Compatible Multipopulation EDO
In ROOT, the main responsibility of the multipopulation EDO is not to find the global optimum but to locate and track multiple moving promising regions. Each subpopulation tracks and covers one promising region. The main purpose of tracking the promising regions for tackling ROOT problems is that the solutions around their summits are more likely to remain acceptable after environmental changes [14], [18], [31]. Using the proposed robustness estimation component's explicit archive, the historical information of the best found positions by each subpopulation in the previous environments is retrieved and used to estimate the robustness of each covered region.
Although many multipopulation EDOs have been developed for tracking the moving global optimum, not all of them are efficient at performing tracking multiple moving promising regions. For example, many state-of-the-art multipopulation EDOs, such as those that form the subpopulations by clustering the individuals based on their fitness and position [32], are defective for this purpose as they may lose track of the inferior (based on fitness) promising regions. Besides, those multipopulation EDOs, whose number of subpopulations and overall population size do not adapt to the discovered promising regions, are not suitable for constructing a ROOT method. Such EDOs are incapable of tracking multiple moving promising regions effectively in the problems whose number of promising regions is larger than the number of the EDO's subpopulations. In such problems, several promising regions cannot be covered due to the limited number of subpopulations.
A suitable class of multipopulation EDOs for performing a stable and reliable tracking of multiple moving promising regions over time are those whose number of subpopulations adapts to the number of discovered promising regions, where the membership of individuals in each subpopulation is fixed and determined based on the individuals' indices [8], [30], [33]. These multipopulation EDOs usually start with one subpopulation (Algorithm 1, line 1) and once it has converged to a promising region, a new subpopulation is initialized (Algorithm 1, line 6). 5 The readers are referred to [2] for a detailed review of the population management and diversity controlling components used in multipopulation EDOs.

B. Proposed Robustness Estimation Component
The first component that is triggered right after an environmental change is the proposed robustness estimation whose main responsibility is to estimate the degree of robustness of the covered promising regions. Estimated robustness values are used for choosing solutions for deployment and also in the population control by the proposed resource allocation. The pseudocode of the proposed robustness estimation component is shown in Algorithm 2.
Each pop i is equipped with an explicit memory M i , which is a circular queue of size m max . After each environmental change, the best found position by pop i in the last environ- . m i is the number of archived solutions in M i . In the case in which the explicit archive is full (i.e., m i = m max ) the oldest archived solution is removed, then g * (t−1) i will be archived. After updating the explicit archive of all subpopulations, a robustness value γ is calculated for each covered promising region (lines 7-9). To this end, for each pop i , γ i is first reset to zero. Thereafter, the acceptability of the archived solutions in M i in the current environment t is evaluated. We first evaluate . This step is repeated until we either reach an archived solution, which is unacceptable or have evaluated all m i archived solutions in M i .
The calculation of γ is costly (i.e., it consumes fitness evaluations), in particular, when the number of discovered promising regions is large. To avoid wasting the computational resources for calculating γ , archived solutions are reevaluated one-by-one from the most recent to the oldest, and once an unacceptable archived solution is detected, the reevaluation process stops. Furthermore, after detecting an unacceptable archived solution, this solution and all older ones are removed from the explicit archive. Actually, when a historical best found solution in a promising region is not accepted in the current environment, it cuts off the chain of robustness in the successive environments. Consequently, we take no account of the older archived solutions in the robustness estimation. We remove the older archived solutions once the chain of robustness is cut off by an unacceptable archived solution (line 12). This pruning mechanism helps to reduce the burden of fitness evaluation consumption caused by using the robustness estimation component.
According to Algorithm 2, γ i ∈ {0, 1, . . . , m i }. The value of γ i indicates the number of successive environments for which if any of the previous best found positions by pop i were deployed, they would have remained acceptable until the current environment. Besides, considering the removal of the unacceptable archived solutions in line 12, all existing archived best found solutions in the tth environment were acceptable from the environment that they were added to the archive until the current environment, i.e., For a pop i , larger values of γ i show that the promising region covered by this subpopulation had some characteristics during the recent γ i environments that resulted in acceptability of the best found solutions over a larger numbers of the environmental changes. Larger values of γ i can be due to several morphological and dynamical characteristics, such as smaller shift severity, smaller fitness fluctuation, and/or wider shape of the promising region. Therefore, larger values of γ i indicate higher likelihood of robustness to environmental changes.

C. Decision-Making Process
A decision maker is responsible to choose a solution for deployment based on γ values calculated by the robustness estimation component (Algorithm 1, line 11). Indeed, γ values, which indicate robustness of the covered promising regions from past to the current environment, form the historical knowledge that is used by the decision maker for choosing solutions for deployment. These solutions are chosen from the promising regions with more reliable and suitable dynamical and morphological characteristics, which are likely robust to the upcoming environmental changes. To identify such promising regions, we use the gathered historical knowledge, i.e., γ values, to estimate their future robustness. We describe the relation between γ values, some dynamical and morphological characteristics of promising regions, and likelihood of future robustness in Section S-I in the supplementary material.
In the proposed ROOT method, the best found position in the promising region with the largest γ value is chosen for deployment. The purpose of choosing such a solution is to maximize the survival time, which is the number of successive environmental changes that the deployed solution can remain acceptable. Applying Algorithm 2 for the promising region covered by pop i , the output indicates that during the last γ i environments, a best found position in this region, which could be successfully deployed, is still reusable until the end of the current environment. Larger γ values for a promising region demonstrate that it has some dynamical and morphological characteristics, which make it more suitable for choosing solutions for deployment that are likely robust to the upcoming environmental changes.
Note that in the proposed robustness estimation component, the archived solutions are used to form the historical knowledge to estimate the future robustness of the promising regions and they are not candidates for deployment. Indeed, similar to [18], the proposed method focuses on robustness of promising regions, which differs from those ROOT methods proposed in [15] and [16], which focus on the robustness of all candidate solutions. This allows us to avoid the complexities and issues of using approximation and prediction components [10], [18], which are necessary for estimating the robustness of candidate solutions used in [15] and [16].

D. Proposed Computational Resource Allocation Component
In most existing multipopulation EDOs, a simple round robin/parallel method is used to allocate computational resources to subpopulations in each iteration [2]. Knowing that subpopulations do not necessarily share the same priority [2], [8], equal allocation of resources to all subpopulations is inefficient [2].
Herein, we propose a new CRA component, which works based on the estimated robustness of the covered promising regions, the role of the subpopulations, subpopulations' progress in their tasks, their age, and the current system status. The proposed resource allocation component uses three thresholds to identify the role of the subpopulations and measuring subpopulations' progress in their tasks. 1) r conv : When the spatial size of a subpopulation is less than r conv , it is assumed that is has converged to a promising region. Otherwise, it has not yet converged to any promising region and is still performing exploration. Such a mechanism is commonly used to identify converged subpopulations in EDOs [2], [34]. In this article, the spatial size λ i of a subpopulation pop i is defined as the Euclidean distance of the farthest pair of individuals [35], which is calculated by 2) r cover : When the spatial size of a subpopulation drops below r cover , we assume that it has converged to the promising region's summit. Note that r cover < r conv and if r cover < λ i < r conv , we assume that although pop i has converged to the promising region, it is still climbing the promising region to reach the summit. 3) r min : This threshold has the smallest value out of the three, i.e., r min < r cover < r conv . When the spatial size of a subpopulation drops below r min , it is assumed that the individuals of pop i have collapsed on the summit of a promising region. In this circumstance, it is assumed that the tracking task has been fulfilled. Many EDOs use r min to deactivate collapsed subpopulations [9], [36]. Algorithm 3 shows the details of the proposed resource allocation component. In each iteration, it determines which subpopulations are allowed to run and use computational resources (Algorithm 1, line 3). As can be seen in Algorithm 3, the proposed resource allocation component is composed of two operational modes, which are selected based on the current solution is no longer acceptable and a new solution must be chosen for deployment. In the following, we describe these two modes.
1) Normal Mode: In the normal mode, the proposed resource allocation prioritizes the following subpopulations. a) Subpopulations with larger γ values: These subpopulations are very important since they are tracking the promising regions that likely contain high quality robust solutions. Consequently, by allocating more computational resources to such subpopulations, their exploitation capability will be accelerated and improved. b) Subpopulations with unfinished exploration/exploitation tasks: The proposed resource allocation component also considers relatively young subpopulations whose γ values are small, but they are important for locating and covering promising regions. A promising region can be considered properly covered when a subpopulation exploits around its summit and tracks it. The last generated subpopulation is the one that is responsible for the vital task of performing exploration in the multipopulation method. This subpopulation fulfills its task once it has converged to a promising region. Besides the explorer subpopulation, the ones that have lately converged to the promising regions, but still have not got close to the summit, are also prioritized by the proposed resource allocation component. The main task of each of these subpopulations is to exploit the promising region and getting close to its summit. Assigning computational resources to such subpopulations to fulfill their current tasks is crucial for providing more accurate information for the robustness estimation component. Besides, the lack of prioritizing such subpopulations may even result in losing the coverage of the newly discovered promising regions.
Allocating computational resources to converged subpopulations is a waste of limited resources. For this reason, resource allocation component stops optimizing subpopulations with lost local diversity. This happens irrespective of the resource allocation component mode, environment number, and γ values. In the resource allocation component, when the spatial size of a subpopulation becomes smaller than r min , it does not get selected until its spatial size is increased by the local diversity control component after each environmental change [2], [10].
In the first environment, the round robin method is used to allocate an identical amount of computational resources to all subpopulations (whose spatial size is larger than r min ) in each iteration since in the first environment, γ = 0 for all subpopulations (Algorithm 3, lines 2-4). After the first environmental change, in each environment where either the previous deployed solution is still acceptable or a new solution has been chosen for deployment, the resource allocation component in normal mode is activated. Below, we describe this process.
In all iterations, the resource allocation component selects the subpopulations with unfinished exploration/exploitation tasks to run in the current iteration. The spatial size of the the last generated subpopulation is always larger than r conv since right after its spatial size falls under the threshold, a new subpopulation is initialized. Therefore, the resource allocation component selects any subpopulation whose spatial size is larger than r conv . Besides, the subpopulations that have recently converged to the basin of attraction of a promising region and are not yet close to the summit are selected. These subpopulations are the ones whose spatial sizes are less than r conv but more than a second threshold r cover . When the spatial size of a subpopulation falls under r cover , we assume that it has converged to the promising region's summit and fulfilled its previous task. Since r cover < r conv , the proposed resource allocation component's normal mode selects all subpopulations with spatial sizes larger than r cover , which are assumed to be not finished with their exploration/exploitation tasks, to run in the current iteration.
The main tasks of the subpopulations with spatial sizes less than r cover include tracking the optimum of the promising region and providing a history of the optimum position in each environment for the robustness estimation component. It is expected that after initial environments and due to prioritizing the subpopulations that have not converged to the summits yet, most existing subpopulations will become tracker subpopulations whose spatial sizes are less than r cover . The reason is that the number of the promising regions in the DOPs considered in this article-and in all existing works-is not excessively high to hinder performance [2]. Consequently, it is expected that after a while, most of the promising regions will be covered by subpopulations that are residing around their summits. Among these subpopulations, the proposed resource allocation component prioritizes those with larger γ values using a probability-based selection process. To this end, at the beginning of each iteration, the resource allocation component assigns a probability p i to each pop i , which satisfies r min < λ i ≤ r cover , i.e., the subpopulations participating in the selection process are the ones whose spatial sizes are larger than r min and smaller than r cover . This probability is calculated as follows: Afterward, a random number is generated in [0, 1] with uniform distribution for each of these pop i and if this number is less than p i , then pop i will be selected by the resource allocation component to execute an internal iteration. In other words, after determining the probability value for each subpopulation using (10), a selection process is independently run for each subpopulation, which decides whether it is chosen to run in the current iteration or not. Since the probability of selecting the subpopulation with the highest robustness value among the participating subpopulations in the selection process is 1.0 using (10), it will definitely be selected. Other participating subpopulations also are selected with probability p i . Thus, at least one subpopulation and at most all participating subpopulations might be chosen. According to this probability-based selection process, the superior subpopulations, i.e., the ones with larger robustness values, have more chance to be selected in each iteration. Note that although the inferior subpopulations are less likely to be selected, they still have a chance to be selected and perform tracking. Moreover, the more frequent the superior subpopulations are selected, the quicker they will be omitted in the selection process due to loss of diversity (once their spatial sizes have shrunk to less than r min ). Consequently, by excluding the subpopulations with larger γ values, the selection probabilities for inferior subpopulations increase.
2) Quick Recovery Mode: When the algorithm responses to an environmental change, if the deployed solution is no longer acceptable, the proposed resource allocation component switches to the quick recovery mode. This mode is designed to ensure that the full potential of the algorithm in finding a better solution for deployment is used, in particular, when the change frequency is high or there is a deadline (temporal constraint) for deploying a new solution [3]. As stated in Section III-C, the best found solution in the promising regions with the highest γ value is chosen for deployment. To this end, we first identify the highest value of robustness among all subpopulations as γ max = max {∀ pop i } (γ i ). Then, we identify the subpopulations whose γ i = γ max and compare their best found solutions. The one with the highest fitness value is chosen for deployment (see Algorithm 4,lines 24 and 25).
In this mode, only the subpopulations whose γ i = γ max run in all iterations while all other subpopulations are hibernated. The resource allocation component remains in this mode until: 1) the spatial sizes of the running subpopulations become less than r min (i.e., sufficient exploitation has been performed);

until stopping criterion is met;
2) the computational budget has run out before the deadline. When one of the aforementioned conditions is met, the resource allocation component will switch back to the normal mode immediately.

E. Detailed Description of an Instantiation of the Proposed ROOT Method
Algorithm 4 shows how the proposed robustness estimation and resource allocation components and those of a multipopulation EDO are assembled to form an instantiation of the proposed ROOT method. Herein, we choose a simple, yet very efficient, multipopulation EDO framework from [9] and [10], which is described in Section S-II in the supplementary material. A complexity analysis of the instantiation of the proposed ROOT method shown in Algorithm 4 is provided in Section S-III in the supplementary material.

IV. EXPERIMENTS AND ANALYSIS
In this section, we first describe the experimental design. Then, we investigate the effectiveness of the proposed robustness estimation and CRA components. We finally compare the performance of a multipopulation EDO equipped with the proposed components and several peer ROOT algorithms.

A. Experimental Design 1) Benchmark Generator:
The experiments in this article are based on various problem instances generated by the generalized MPB (GMPB) [28], [37]. GMPB is a benchmark generator, which is capable of generating landscapes with a controllable number of promising regions with a variety of parametric and changeable characteristics, including symmetry, condition number, irregularity, roughness, modality, and variable interaction. GMPB has several parameters (see Table I) that can be set by the user to generate problem instances with a vast variety of morphological and dynamical characteristics and difficulty levels. Detailed information of the GMPB used in this article is provided in Section S-IV in the supplementary material. The MATLAB source code of this benchmark is available from [38].
2) Performance Indicator: The focus of this article is on solving ROOT problems in which the main goal is to find solutions for deployment in order to maximize the average number of environments that the deployed solutions remain acceptable. It is worth to mention that none of the algorithms examined in this article, including the proposed one, are designed to maximize the fitness of deployed solutions or minimize the switching cost, which are related to other classes of ROOT problems or can be considered as other objectives [11], [12], [17]. Considering the focus of this article, we compare the performance of the algorithms according to the average survival time [17], which is the most commonly used performance indicator in the field [16], [18].

3) Algorithms:
In the experiments, we use the multipopulation framework from [10] for all multipopulation-based methods, including the proposed one. We also use PSO with a constriction factor [39] as the optimization component in the multipopulation framework. To evaluate the effectiveness of proposed resource allocation and robustness estimation components, we add them to mPSO. This approach is denoted as mPSO +CRA +RE in the experiments. Besides mPSO +CRA +RE , we also use three other multipopulation comparison algorithms, which are: mPSO, mPSO b , and mPSO sh . mPSO is a tracking the moving global optimum method that is adapted to tackle ROOTs in which the best found solution (in terms of fitness) is chosen for the next solution for deployment. mPSO Rb uses the ROOT decision maker from [20], where the best found solution among the candidate reliable promising regions is chosen for deployment. mPSO Rsh applies (7) (for which estimated shift and height severity values are needed) to choose the next solution for deployment [18].
We also choose the fitness prediction-based ROOT methods from [17] (denoted as PbM F ) and [15] (denoted as PbM J ) as comparison algorithms. For the experiments, we assume that these methods have access to the previous environmental parameters; thus, they do not need any approximation method. It should be noted that the fitness evaluations used for training the predictor in the previous environments are not counted toward the overall computational cost. This means that the results obtained by these algorithms are not affected by the approximation error and can therefore be taken as an upper  [38] bound for what the algorithms are capable of achieving in practice.
In this article, we focus on DOPs with visible environmental changes where the optimization algorithms are informed about environmental changes by other parts of the system, such as sensors and agents, similar to many real-world DOPs [2]. Therefore, the examined algorithms in this article do not use any change detection component. If needed, a simple reevaluation-based change detection component [40] could be added to the algorithms. Table I shows the parameter settings of GMPB. The experiments are done on the problem instances with various numbers of promising regions (m), acceptability threshold (μ), and dimensions (d). The parameter settings chosen in Table I are commonly used in the ROOT and EDO literature. Different parameter settings result in problem instances with different difficulty levels. In addition, different parameter settings generate problem instances with different maximal robustness values for both problems and promising regions. By increasing μ, the regions containing the robust solutions shrink, the maximal possible survival time decreases, and finding robust solutions becomes more challenging [16]. The higher the number of promising regions m, the easier to find robust solutions. This is due to the fact that in landscapes with larger numbers of the promising regions, the size/number of areas containing robust solutions with better qualities increases [18]. In addition, in problems with larger m values, the promising regions are likely to overlap and support the deployed solution to remain acceptable for further environments. Finally, by increasing the dimension, the problem instances become more challenging.

4) Parameter Settings:
We use the parameter settings extracted from the sensitivity analysis in [10] for our multipopulation EDO framework and also for PSO in the multipopulation methods. Besides, the parameter settings of the proposed components are taken  Table II. For the parameter settings of the comparison algorithms, the values suggested in their original references are used. Our investigations also indicate that these algorithms show their best efficiency with those suggested settings.

B. Experimental Results
The statistical results provided in this section are based on 31 independent runs. For statistical analysis, we use multiple comparison tests using the Wilcoxon rank-sum test with Holm-Bonferroni p-value correction and α = 0.05.

1) Effect of the Proposed Components on Performance:
In this section, we investigate the effectiveness of the proposed robustness estimation and CRA components. To this end, we compare the performance of mPSO, mPSO with robustness estimation component (mPSO +RE ), and mPSO with both proposed robustness estimation and CRA components (mPSO +CRA +RE ). Fig. 1(a) compares the average survival time over time by these three methods on GMPB with the default parameter settings from Table I. We also compare the average percentage of the previously deployed solutions that are reused (i.e., remaining acceptable) by these methods in each environment in Fig. 1(b). Comparing the performance of mPSO and mPSO +RE in these plots based on average survival time and the average percentage of acceptability of the previously deployed solutions demonstrates the effectiveness of this component in choosing more robust solutions for deployment.
As can be seen in Figs. 1(a) and (b), adding the proposed CRA component further improves the performance of the algorithm. This improvement is a result of systematic control of the computational resource consumption by each subpopulation according to their estimated robustness value (provided by the proposed robustness estimation component), role, task achievement, and current system status. To further investigate the effectiveness of the proposed resource allocation, we compare the performance of mPSO +RE and mPSO +CRA +RE on the problem instances generated by GMPB with different numbers of the promising regions and the default settings for the remaining parameters. The results are compared in Fig. 2. As can be seen in this figure, by adding the proposed resource allocation component, the performance is improved, in particular, in the problem instances with larger numbers of the promising regions. In such problems, larger numbers of subpopulations are generated to cover the promising regions. Thus, the role of resource allocation becomes more vital since larger numbers of the fitness evaluations are needed in each iteration. mPSO +RE suffers from a shortage of the available computational resources before the deployment time. However, for the mPSO +CRA +RE , this challenge is ameliorated by systematic allocation of the computational resources to the subpopulations with higher robustness and taking the system status into account.
As described in Section III-D, besides prioritizing the subpopulations with larger γ values, the proposed CRA component also prioritizes the subpopulations, which are performing exploration or recently have converged to a promising region and are moving toward the summit. This systematic Fig. 2. Investigating the effect of the proposed CRA by comparing the obtained average survival time by mPSO +RE and mPSO +CRA +RE in problem instances generated by GMPB with different numbers of promising regions (m) and the default parameter settings from Table I for the rest of the parameters. Fig. 3. Investigating the effect of the proposed CRA by comparing the ability of finding and covering/tracking promising regions in mPSO +RE and mPSO +CRA +RE over 100 environments of GMPB with 50 promising regions (m = 50) and the default parameter settings from Table I for the rest of the parameters. The plots are obtained by averaging the results of 31 independent runs. Note that usually some smaller promising regions are covered by larger ones, thus, the number of visible promising regions is usually less than m. Besides, by changing the size and location of promising regions, the number of visible promising regions changes over time [28].
prioritizing approach considerably improves the abilities of exploration, exploitation, and tracking, which results in discovering and covering larger numbers of promising regions. Fig. 3 shows the effect of using the proposed resource allocation in the ability of the algorithm in finding and covering/tracking promising regions over time. Note that a promising region is considered covered if there is a subpopulation whose individuals reside in the basin of attraction of the promising region. As can be seen in this plot, thanks to the proposed resource allocation component, mPSO +CRA +RE covers larger numbers of promising regions in comparison to mPSO +RE . By covering larger numbers of promising regions in mPSO +CRA +RE , the possibility of missing promising regions, which may contain solutions with higher robustness, decreases that in turn results in improving the performance of the algorithm in finding better robust solutions.
2) Comparison With Peer Algorithms: In this section, we compare the results obtained by mPSO +CRA +RE and the peer algorithms described in Section IV-A3 in solving problem instances generated by GMPB with different dimensions d, numbers of promising regions m, and acceptability threshold values μ. The average survival time (and standard error) obtained by the algorithms is reported in Table III, where  does not suffer from inaccuracy in estimating these values, especially, in the DOPs with random and/or heterogeneous dynamics [2] where estimating these parameter values is error prone. Besides, additional morphological and dynamical characteristics of the promising regions are also implicitly considered in the calculation of the estimated robustness values in mPSO +CRA +RE , which are not taken into consideration for choosing solutions in mPSO Rsh and mPSO Rb .
The proposed CRA component is another major feature of mPSO +CRA +RE . The used resource allocation methods in mPSO Rsh , mPSO Rb , and mPSO are designed for tracking the moving global optimum [10], and they do not consider any attribute related to the robustness of promising regions in controlling subpopulations. On the other hand, the proposed resource allocation is particularly designed for tackling ROOT problems, which controls the subpopulations based on several factors, including the robustness of the promising regions and the acceptability of the deployed solution. Besides, the proposed resource allocation takes the roles of subpopulations, their convergence status, and their task achievements into consideration.
To further investigate the performance of mPSO +CRA +RE in problem instances with different characteristics, we carry out additional experiments on the problem instances generated by GMPB with different acceptability ROOT threshold values, change frequencies, shift severity degrees, and computational budget values. The results are reported in Section S-VI in the supplementary material. Moreover, to show the independence of the proposed components with respect to the optimization component used, we compare the performance of the multipopulation methods when they use differential evolution (DE) [42] as the optimization component in Section S-VII in the supplementary material. The results indicate the superiority of the proposed ROOT method when DE is used as the optimization component.

V. CONCLUSION
In this article, we have presented two new componentsrobustness estimation and dual-mode CRA-for ROOT methods. The robustness estimation component is responsible for estimating the robustness degree of the covered promising regions, while the systematic dual-mode robustness-based CRA component controls the subpopulations. These two components and those of a multipopulation EDO, which is capable of tracking multiple moving promising regions, are assembled to form a new ROOT method. Unlike the existing components for determining robust solutions, the performance of the proposed robustness estimation component does not rely on the oversimple characteristics of the benchmark problems, such as fixed relocation length of the promising regions, unimodality and smoothness of the promising regions, and/or low dimensionality of the problem. Moreover, using the proposed dual-mode CRA component, the proposed ROOT method is the first one that takes robustness and the system status into account for controlling the subpopulations. The proposed ROOT method and a set of peer methods have been used for maximizing the average survival time of the deployed solutions in 48 different problem instances with various characteristics generated by GMPB. The experimental results have shown the superiority of the proposed method in almost all test cases.
In our algorithm, we monitor the spatial size of subpopulations to determine their role and also progress in carrying out their tasks, where we used a simple widely used method [35] for calculating the spatial size of subpopulations. Designing a more advanced and systematic spatial size monitoring method that takes some problem characteristics, such as illconditioning and asymmetry into account, is a potential future research direction.
In the proposed robustness estimation component, the value of m max should be set considering the scope of robustness in the problem. In this article, the value of m max is fixed over time, which may not be efficient for solving heterogeneous DOPs [2] in which dynamical behavior changes over time. In such problems, the scope of robustness significantly changes over time following the changes in the dynamical characteristics, such as change frequency and severity. A potential future work will be designing parameter adaptation mechanisms [43] for adapting the value of m max to the changing scope of robustness over time in heterogeneous DOPs.
There are many real-world problems whose search spaces contain multiple moving promising regions and the multipopulation-based methods, such as our proposed method, are efficient in solving them. A potential future work will be solving a real-world ROOT problem. An example of real-world ROOT problems is crowd monitoring and management [44], which is a dynamic covering location problem [45]. In this problem, the locations of security agent units are changed over time based on the status of the crowd. However, frequently changing the locations of security agent units is undesirable as it disturbs the monitoring task. Consequently, in this problem, we seek solutions (i.e., the locations of the security agent units) that can remain acceptable for a longer time.