Pareto Set Prediction Assisted Bilevel Multi-objective Optimization
Bing Wang Hemant K. Singh Tapabrata Ray












bing.wang@unsw.edu.au

  School of Engineering and Technology, University of New South Wales Canberra


 Australia




h.singh@unsw.edu.au

  School of Engineering and Technology, University of New South Wales Canberra


Australia




t.ray@unsw.edu.au

  School of Engineering and Technology, University of New South Wales Canberra


Australia















§ ABSTRACT 
 Bilevel optimization problems comprise an upper level optimization task that contains a lower level optimization task as a constraint. While there is a significant and growing literature devoted to solving bilevel problems with single objective at both levels using evolutionary computation, there is relatively scarce work done to address problems with multiple objectives (BLMOP) at both levels. For black-box BLMOPs, the existing evolutionary techniques typically utilize nested search, which in its native form consumes large number of function evaluations. In this work, we propose to reduce this expense by predicting the lower level Pareto set for a candidate upper level solution directly, instead of conducting an optimization from scratch. Such a prediction is significantly challenging for BLMOPs as it involves one-to-many mapping scenario. We resolve this bottleneck by supplementing the dataset using a helper variable and construct a neural network, which can then be trained to map the variables in a meaningful manner. Then, we embed this initialization within a bilevel optimization framework, termed Pareto set prediction assisted evolutionary bilevel multi-objective optimization (PSP-BLEMO). 
Systematic experiments with existing state-of-the-art methods are presented to demonstrate its benefit. The experiments show that the proposed approach is competitive across a range of problems, including both deceptive and non-deceptive problems [This is author submitted version of the work currently undergoing peer-review].







<ccs2012>
<concept>
<concept_id>10002950.10003714.10003716.10011136.10011797.10011799</concept_id>
<concept_desc>Mathematics of computing Evolutionary algorithms</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>
<concept_id>10003752.10003809.10003716.10011138</concept_id>
<concept_desc>Theory of computation Continuous optimization</concept_desc>
<concept_significance>500</concept_significance>
</concept>
</ccs2012>

[500]Theory of computation Continuous optimization
[500]Mathematics of computing Evolutionary algorithms













Pareto Set Prediction Assisted Bilevel Multi-objective Optimization 
    Tapabrata Ray
    
====================================================================





§ INTRODUCTION

A number of real-world problems manifest as a hierarchical optimization problem, where the objective(s) of an upper-level (UL) problem are optimized subject to the optimality of a lower-level (LL) problem. The UL and LL are also referred to as leader and follower tasks, respectively, and may contain their own additional constraints. Such problems are referred to as bilevel programming or bilevel optimization problems (BLP) <cit.>; and are the most-studied subset of a more generalized class of multi-level optimization problems <cit.>. Such problems are of interest in a number of application domains <cit.>, including transportation<cit.>, economics <cit.>, engineering <cit.> and management <cit.>, to name a few. Moreover, they also pose certain unique theoretical and algorithmic challenges <cit.> compared to traditional optimization (referred to occasionally in this study as single-level optimization for specificity). These include, for example, being NP-hard even when the problems are linear at both levels <cit.>, deceptive nature of UL evaluations for solutions that are sub-optimal at LL, and challenges in environmental selection and benchmarking strategies <cit.>. Considering the above, bilevel problems have attracted interest from both academic researchers and practitioners. Formally, bilevel problems can be represented as shown in Eq. <ref>. 


    [ 𝐱_u∈𝕏_uMinimize F(𝐱_u,𝐱_l) = F^1(𝐱_u, 𝐱_l), ...., F^M(𝐱_u, 𝐱_l);                                                     Subject to:;     𝐱_l ∈𝐱_l∈𝕏_largmin f(𝐱_u, 𝐱_l) =f^1(𝐱_u,𝐱_l) … f^N(𝐱_u,𝐱_l);     subject to:  g^1(𝐱_u, 𝐱_l),g^2(𝐱_u, 𝐱_l)… g^q(𝐱_u, 𝐱_l) ≤ 0;                 G^1(𝐱_u, 𝐱_l), G^2(𝐱_u, 𝐱_l)… G^p(𝐱_u, 𝐱_l) ≤ 0;                                                                 ]

Here F^1 … F^M and f^1 … f^N are objective functions at UL and LL, respectively. The UL and LL design variables are denoted using {𝐱_u, 𝐱_l}, sampled from their design spaces {𝕏_u, 𝕏_l}, respectively. The function g(𝐱_u, 𝐱_l) represents a constraint for the LL problem, while G(𝐱_u, 𝐱_l) determines feasible space for the UL problem. For non-exact techniques, the equality constraints, if they exist, are often converted into inequalities with a tolerance, hence equality constraints are not included in Eq. <ref> for brevity. 

Bilevel optimization problems have their origins in Stackelberg game theory and economics <cit.>, but have been since studied in several other fields. However, overwhelming majority of these studies have been dedicated to problems where both levels have a single objective (BLSOP), i.e., M=N=1. A number of competitive algorithms have been developed for this class of problems utilizing ideas such as problem transformation <cit.>, hybridization <cit.>, surrogate-assisted search <cit.>, transfer learning <cit.>. For more comprehensive overview of the works in the area of BLSOP, the readers may refer to review papers such as <cit.>. 

However, there has been relatively scarce attention paid to bilevel problems with multiple objectives at one or both levels. The bilevel multi-objective optimization problems (BLMOPs), also denoted with acronym MOBO in <cit.>, typically have multiple objectives at UL, while the LL may have one or multiple objectives. The special case where the UL has a single objective, but LL has multiple objectives, is referred to as semi-vectorial bilevel optimization (SVBO) problem <cit.>. A number of real-world problems can be formulated as BLMOPs, and thus could benefit from efficient methods to solve such problems. In <cit.>, an environmental economics problem is discussed wherein a company wants to set up a gold mine in a specific geographical region. The government acts as the UL decision-maker in this case, with the objectives of maximizing revenues generated by the project (tax, jobs etc.) and minimizing the detrimental impacts on the environment. The mining company, as the LL decision-maker, aims to maximize its profits as well as maximize its reputation/public image. Another problem discussed in <cit.> is that of hierarchical decision-making in a company, wherein the objectives of a chief executive officer (CEO) are to maximize the quality of the products and the company profits. On the other hand, the branch heads working at a lower level seek to maximize the branch profits and the worker satisfaction. In the context of border security, the problem discussed in <cit.> involves UL objectives of maximizing total weighted exposure of the minimal exposure path, minimizing sensor relocation time and the number of sensors relocated, while the LL seeks to minimize the total expected weighted exposure of the intruder’s minimal exposure path. Some other BLMOP application examples include  transportation planning and management <cit.>, manufacturing <cit.>, machine learning problems <cit.>. A number of additional applications from the domains of logistics, environmental economics and manufacturing can also be found in <cit.>.

In this study, we are predominantly concerned with addressing BLMOPs where both levels have multiple objectives. Given the limited work in the field so far, the scope of study is limited to two objectives, but in principle higher number of objectives fall within the same category. Beyond the challenges already encountered for BLSOPs, BLMOPs have their own characteristic challenges. The key challenge is that for any given UL solution, the optimum solution at LL is not just a single solution but a set of trade-off solutions (Pareto front, PF). This situation occurs rarely in BLSOPs, where the LL optimum has multiple global optimum solutions, but is ubiquitous for BLMOPs. This implies that for each 𝐱_u, not only a significant number of function evaluations are required at LL to find a good PF approximation, a large set of solutions (corresponding to LL PF) also needs to be evaluated at UL. This results in proliferation of function evaluations used by approaches, especially those that rely on nested strategies. The second, related challenge is that of autonomy in decision-making by UL and LL. Since LL decision-maker has an entire PF approximation to chose from for sending to UL decision-maker, co-operation or conflict between the two decision-makers may result in sub-optimal PF at the UL. In most of the existing works, an optimistic version of the problem is assumed, wherein the UL decision-maker has the authority to choose solutions from the LL PF that bring UL the most benefits. In what follows, we highlight some of the representative approaches that address BLMOPs. Moreover, some of the works convert the LL problem from multi-objective to single-objective through the use of value functions in order to provide a single solution to the UL <cit.>. For a more extensive coverage of literature in BLMOPs, the readers are referred to  <cit.>.

Classical/analytical techniques:  
Much like the case of BLSOPs, some of the early efforts were directed towards solving BLMOPs utilizing mathematical techniques for exact solutions. The pre-condition to the application of these techniques is often that the underlying response functions (objectives, constraints) need to satisfy certain regularity conditions. For instance, in <cit.>, an approach is developed to solve BLMOPs where all objectives and constraints at both levels are linear. Three different cases were solved with regards to anticipated behavior of the LL decision-maker, namely, optimistic, pessimistic, and historical. In <cit.>, bi/multi-level problems with multiple objectives are solved under a scenario where the decision-makers have certain tolerances, modeled using fuzzy set theory. In <cit.>, an interactive method was presented for BLMOP by replacing the LL problems with Kuhn–Tucker conditions. In <cit.>, an approach to solve non-linear BLMOPs was presented. Differentiability is generally assumed for calculating gradients in such approaches.
In <cit.>, BLMOP is linearly scalarized to generate Pareto frontier solutions, and then a filtering mechanism is constructed to maintain distribution for these solutions. Linear functions are assumed for the given problems.  In <cit.>, Pareto solutions of bilevel linear problems have been characterized after converting bilevel problem into single level mixed 0-1 programming problems. 
Moreover, in most of these works, the working principles are demonstrated numerically on a very limited set of problems (often one or two), with small number of variables. 

Metaheuristic techniques:  In order to handle cases where the underlying functions do not satisfy the required mathematical properties, or are not available altogether (so-called “black-box” functions), metaheuristic techniques, such as evolutionary algorithms or swarm intelligence algorithms, are often chosen. However, this versatility typically comes with the requirement to evaluate large number of candidate solutions to reach near-optimal solutions. Given the structure of a bilevel problem, a common way to solve it using metaheuristic methods is through a nested approach, where a standard multi-objective evolutionary algorithm (MOEA) or equivalent methods are used at both levels. Towards this end, non-dominated sorting algorithm II (NSGA-II) with a special population structure was utilized in <cit.> in a nested manner to solve BLMOPs. This was further hybridized with local search methods in <cit.>. Differential evolution was implemented in a nested form for solving BLMOPs in <cit.>.
With similar intent, nested search has also been proposed using other metaheuristic approaches, such as particle swarm optimization (PSO) and simulated annealing  (SA) <cit.>.
Recently, a variable decomposition based cooperative co-evolutionary method was proposed to improve search efficiency of nested search structure of BLMOP <cit.>. 
A generic framework to solve BLMOPs was also presented recently in <cit.>, wherein some gain in efficiency was realized via a representation that utilizes grouping of certain variables as families. 

Since the implementation of meta-heuristic search in a nested mode involves, in particular, a large number of LL function evaluations (FE), a pertinent effort is to develop strategies to expedite the search at LL and reduce LL FE. Towards this end, transfer of solutions obtained corresponding to a neighboring UL solution as a seeding population has been investigated for BLSOPs in <cit.>,<cit.>,<cit.>, and more recently, extended for BLMOPs in <cit.>. Though the results were promising, the performance of such simple transfer strategies is subject to neighboring landscapes being similar, as was demonstrated for BLSOPs in <cit.>. Another promising direction to improve LL search efficiency is by generating LL optimal solutions through certain prediction mechanism. Often, the LL optimal solutions are not entirely independent of each other, since they are related by the UL variables that act as parameters for the LL problems. Consequently, by learning from the dataset of previously evaluated solutions, predictions can be made regarding the optimal 𝐱_l for a new candidate 𝐱_u. The predicted LL optimum can be used as a starting point for local search, or to seed initial population to reduce the effort required in carrying out the LL optimization. In mf-BLEAQ <cit.>, a quadratic fiber based surrogate model was proposed to predict LL Pareto sets (PS). Using predicted solutions to start the LL search provided a head-start. However, this modeling scheme needs to build a model for each variable of a solution, and for LL PS there is a need for multiple solutions to cover PF. Therefore, a large number of models are required to be built and maintained. A more recent study <cit.> proposes to reduce the number of prediction models for LL PS to a single model, i.e., conditional generative neural network (cGAN). Relying on the approximation capacity of cGAN, the proposed algorithm cG-BLEMO only needs to maintain one model. However, due to Gaussian noise introduced in cGAN, LL solutions generated by cGAN do not follow the distribution of the PS/PF closely, requiring additional search to improve the solutions to capture LL PS. Furthermore, to keep the training time practicable for cGAN, the datasize used had to be kept relatively small in the implementation (typically up to 800).

Continuing the above line of inquiry, in this study, we aim to advance the research direction that uses LL PS prediction for solving BLMOPs efficiently. Towards this goal, we build a simple yet accurate model customized for BLMOPs, that can be used to predict LL PS in lieu of evolutionary search to drive majority of the search. The key contributions of this work can be summarized as:

    
  * Ordinarily, the LL optimum comprises multiple solutions (PS) for each UL solution 𝐱_u. Thus, the dataset of previously evaluated solutions results in a one-to-many mapping with regards to LL PS, which is unsuitable for building prediction models. We propose a simple transformation of the data to make them amenable for creating one-to-one mapping. 
    
  * We then use simple feedforward neural networks to build the prediction model that maps 𝐱_u and optimal 𝐱_l (LL PS). The model training is relatively fast, uses fewer hyperparameters compared to more sophisticated models previously used in the literature, and is more accurate in terms of its LL PS prediction. 
    
  * We embed the above PS prediction into the nested search for BLMOPs. After accumulation of training data and construction of the NN model, LL search can gradually be delegated to prediction. To maintain/improve the prediction accuracy, NN is re-trained periodically during the run with new evaluated solutions.
    
  * We benchmark the proposed approaches on a diverse set of BLMOPs from the literature, including deceptive problems that have been scarcely studied. Comparison with state-of-art algorithms highlights its competent performance in terms of proximity to optimum, and efficiency in terms of function evaluations used.   

 

The remainder of this paper is structured as follows. In Section <ref>, we present the main proposed idea for PS prediction and show proof-of-concept results to establish its viability. 
Then, this prediction method is embedded within a nested bilevel multi-objective optimization framework, the details of which are presented in Section <ref>. The numerical experiments and discussions are presented in Section <ref>, followed by concluding remarks in Section <ref>. 




§ THE MAIN IDEA AND PROOF-OF-CONCEPT




In standard (single-level) multi-objective optimization, it is common to seek intermittent estimation of Pareto front/Pareto sets to aid the search, especially when the underlying response functions (objectives, constraints) are computationally expensive in nature. This can be done by training a surrogate model on an archive of evaluated solutions to establish a mapping between the input 𝐱_i and output response(s) y_i, as shown on the left side in Fig. <ref>. Thereafter, an internal optimization can be conducted using the surrogate model to come up with solutions that are likely to be Pareto optimal, which can then undergo the expensive evaluation <cit.>. Given that the accuracy of the models is dependent on the dataset, this process usually needs to be conducted multiple times over a run. More recently, there have also been attempts <cit.> to predict the Pareto optimal front based on preference vectors. This is done by utilizing preferences (quantified using, e.g., reference vector coefficients) with regards to different objectives as inputs, and certain known Pareto optimal solutions as output. The mapping between these can be done using various machine learning models such as neural networks, Gaussian process or hypernetworks. The mapping can then be used to generate new solutions that are likely to be Pareto optimal; which can be used, for example to increase the density of solutions on the PF, or aid the search interactively.

For the above modeling tasks, the dataset usually has a one-to-one or many-to-one relationship between inputs and outputs, for which there is no mathematical ambiguity in terms of building the mappings. In BLMOPs however, this is mostly not the case. If the models are built individually for LL or UL functions to aid the search independently, the above models and types of mappings can still be maintained. However, to potentially bypass the LL optimization for some UL candidate solutions entirely by predicting the LL PS, a different type of mapping is required. As shown in the right side of Fig. <ref>, the available dataset in this case may be the previously evaluated solutions, 𝐱_u^i. For each of these solutions, a lower level PS approximation has been found using an MOEA, denoted as {𝐱_lj^i*,j=1… m }. That task is to predict the PS {𝐱_l^* } directly for a new 𝐱_u sampled during the UL search. For BLSOPs, in most cases this can still be  modeled as a one-to-one mapping <cit.>, but for BLMOPs, it is evident that such a mapping is not straightforward due to the same value of input 𝐱_𝐮 for multiple values of 𝐱^*_l. 






To address some of the issues above, we propose a simple method in this work to transform the data and use neural networks (NN) for creating a mapping between 𝐱_u and 𝐱^*_l. Neural networks can inherently support vectors as outputs (which is the case with 𝐱^*_l), and are also proven to be universal approximators given sufficient training data. Given the number of points accumulated with LL PS corresponding to multiple 𝐱_u values, the size of the dataset works in favor of neural networks for BLMOPs to create accurate mappings. Of course, the issue of one-to-many mapping still needs to be resolved. We use a helper variable to supplement the input data to overcome this. Our proposed approach is quite simple in its structure and discussed below with an example. 

To begin with, the dataset that will be provided to the NN needs to be generated and conditioned. To generate the initial data, the LL PS needs to be identified for certain initial samples of 𝐱_u. These LL PS (or approximations thereof) can be obtained using an MOEA at the LL. Once obtained, the available data resembles the one shown on the right in Fig. <ref>. Then, the data for each 𝐱_u is conditioned schematically as shown in Fig. <ref>. To do so, the LL PS approximations {𝐱_l^*} are first sorted based on increasing order of one of the LL objectives; chosen (without loss of generalization) as f^1 in this case. This ordered set forms the output values for the NN. For the inputs, the value of a given 𝐱_u is replicated m times, where m is the number of solutions in the LL PS approximation corresponding to 𝐱_u. In addition, a helper variable r is appended to the inputs. The variable r assumes a uniformly spaced values in the range of 0 to 1, where r_1 corresponds to 0 and r_m corresponds to 1. The r values can thus be generated easily as per the number of points in the LL PS approximation. As it can be seen, the resulting dataset has unique input vectors, each mapping to a single output vector 𝐱^*_l. Therefore, the dataset is suitable for many-to-one or one-to-one mapping. The significance of having 𝐱^*_l in the sorted order of f^1, which corresponds to sorted order of r in the inputs, can be inferred as providing a sense of direction. This can be thought of as akin to mapping solutions along different preferences or reference vectors in, e.g., decomposition-based evolutionary algorithms. It is worth noting that, using a single helper variable r is only suitable for solving bi-objective problems, because there is an unambiguous spatial ordering among the PF points for bi-objective problems. That is, the PF points can be traversed from one corner of the PF to another, corresponding to the ordered set of reference vectors based on the monotonic progression of r values. Such ordering is not straightforward for problems with more than two objectives. Using the proposed approach for more than two objectives may be feasible following similar ideas of decomposition, but will require (a) more than one helper variables, and (b) careful consideration of the ordering of the data points and helper variables. Given that most of the existing BLMOPs involve two objectives <cit.> and that the existing methods require significant developments even for bi-objective problems, we limit the scope of this study to two objectives.


Once the above processing has been applied to all available 𝐱_u and the corresponding 𝐱^*_l, the data is stacked to form the dataset for training a feed-forward NN (FNN). It takes the inputs (𝐱_u,r) and builds a model 𝕄' to predict the output 𝐱^*_l. This model can then be integrated in the bilevel search method to predict the LL PS for any new candidate 𝐱_u under consideration, or can form a seed population to expedite the LL search towards its PS. 

The benchmark problem DS2 <cit.>, with two variables at each level, is used here for proof-of-concept. The theoretical LL PS can be inferred for the problem for any given 𝐱_u, which can be used to verify the predication accuracy. For this illustration, we sample 10 random solutions (𝐱_u) at the UL. For each 𝐱_u, we sample 20 Pareto solutions (𝐱^*_l) on the LL PS. This data is then supplemented with the helper variable values as discussed above, and the FNN is trained. 

We use MATLAB feedforward neural network tool () to create the mapping from {(𝐱_u, r)} to {𝐱^*_l}. One hidden layer with 4 nodes is set for this problem. Default parameter settings are used for the training step, mean square error is used as loss function, and backpropagation is used for training (More detailed parameters will be discussed in Section IV). Next, we generate a new 𝐱_u and attempt to predict its PS. The generated UL variable, 𝐱_u=[1.2, 1.2] is combined with  r values  to predict 20 LL PS solutions. The image of predicted PS in the LL objective space is shown in Fig. <ref>. 



It can be seen that predicted solutions lie on the true PF as well as exhibit a good diversity on it. In Fig. <ref>, in variable space, we plot PS and predicted solutions. It shows consistency between PS and predicted solutions, For comparison, we also show the results obtained when the r is not appended to 𝐱_u in a uniformly increasing sequence but randomly generated instead (this  is also equivalent to not sorting the 𝐱^*_l based on f^1). The results can be seen in  Fig. <ref> and Fig. <ref>, where the mapping is able to predict solutions close to only a small part of the PF and PS, respectively. Lastly, we present 20 solutions that are randomly sampled in the LL search space for the same 𝐱_u. This is shown in Fig. <ref>, and highlights that if the LL search was conducted in a standard manner starting from a random population, the initial solutions will be far away from the PF. Instead, if the predicted PS such as that in Fig. <ref> is utilized, the search can be bypassed entirely for some generations, or seeded on/close to the PS, expediting the LL search and potentially saving significant number of LL evaluations. Using this central idea, we build an algorithmic framework to solve BLMOPs, termed Pareto set prediction assisted bilevel evolutionary multiple objective optimization algorithm (PSP-BLEMO) next. 




§ PROPOSED ALGORITHM



The general framework of PSP-BLEMO[The code and data will be made available for research purpose after the review process.] is outlined in Algo. <ref> and Fig. <ref>, followed by description of its key components in the following subsections. At UL, the search is conducted through an MOEA, while the LL search has the PSP model integrated in it. 










 §.§ UL search and variable association check

The UL search commences with an initial LL variable association check, followed by the creation of an initial population P_1
generated through a random uniform distribution. An LL variable association check is done for accommodating an uncommon feature of some BLMOPs, referred to as variable association ambiguity (VAA) <cit.>. If a BLMOP has LL VAA, a subset of its LL variables participate solely in UL objectives. A detailed discussion of VAA can be found in <cit.>, here we only briefly introduce the VAA checking step adapted to this study, to keep the primary focus of this study on PSP module. The VAA check is done by perturbing an LL variable while holding all other variables constant. By comparing the LL objective values before and after the perturbation, it becomes possible to ascertain whether this variable is associated with LL objectives. Likewise, the UL objective values are also checked before and after the perturbation to ascertain which LL variables affect them. This process is repeated for all LL variables one by one. Subsequently, a binary vector v (Algo <ref>, Line <ref>) is generated to record the association status of each variable: 0 signifies no VAA, while 1 indicates otherwise.[Note that theoretically another case is possible where an LL variable is completely redundant, i.e, does not feature in LL or UL. We have omitted discussion of such cases for brevity in this study, but they can be easily filtered out from the search entirely through the above VAA checks.]
Upon completion of the VAA check and the generation of P_1, we proceed to evaluate P_1. To evaluate each individual 𝐱_u^i within P_1, LL search is conducted using an MOEA to identify the corresponding LL PS approximation {𝐱_l^i*} (Algo. <ref> Line <ref>). In case {𝐱_l^i*} is empty, such a solution is assigned a ∞ value at UL. The special cases of VAA require additional search components, which will be discussed later in Algos. <ref>-<ref>.




The 𝐱_u of the non-dominated (ND) solutions in P_1 and the {𝐱_l^*} corresponding to each of them are paired and stored in an archive (Algo. <ref> Line <ref>). All solutions that went through LL search are accumulated into the archive for training/updating the PSP model (Algo. <ref> Line <ref>). 
As indicated in the previous section, the helper variable r is also generated and appended to 𝐱_u in the archive. After the first population P_1 has been evaluated through LL search, UL offspring population is generated by using evolutionary operators (Algo. <ref>, Line <ref>). Since the PSP model has been trained in the previous step, it can be used to supplement LL search for evaluation of the UL new offspring. This can be done by initializing the LL population using the predicted PS (Algo. <ref>) for the LL search. Here, we introduce an NN update parameter γ to explore the potential of the predictor, and  a parameter ds for training data size threshold. Without sufficient training data, the accuracy of the model may not be good enough to reliably skip LL search. Therefore, LL search is conducted as usual, until training data size meets ds (Algo. <ref>, Line <ref>). Thereafter, the LL search is only run every γ generations, initialized with PSP. The LL search is skipped in the intermediate generations entirely (Algo. <ref>, Line <ref>). In this way, the LL search can be forgone for most of the generations to save on evaluations. The training data for NN are only accumulated from solutions for which LL search is conducted to ensure the accuracy of the evaluations in the dataset (Algo. <ref>, Line <ref>). When the number of data points in the archive exceeds ds, the most recent ds data points are used for building the model, to keep the model more accurate in the current region of search.
After each UL child has received its corresponding LL solutions {𝐜_l^i*} and evaluated at UL (Algo. <ref>, Line <ref>), UL ND archive is first updated. Then parent and child populations are combined and sorted according to feasibility first (FF), non-dominance (ND) and crowding distance (CD) (Algo. <ref>, Line <ref>). The next generation UL population is selected from combined sorted population using environmental selection mechanism outlined in Algo. <ref>. This evolution continues until termination criterion is met, and the archive of ND solutions obtained is returned. 




 §.§ PSP and PSP-assisted LL search

The role of PSP is to make LL search efficient by either replacing it altogether with predicted PS, or providing LL search a better starting population than random initialization. For LL PS prediction, when a new 𝐱_u is given to PSP, a list of linearly spaced helper value r (r∈[0, 1]) is first generated (Algo. <ref> Line <ref>). Then 𝐱_u is replicated to achieve the same number (n^l) of solutions that exist in r (Algo. <ref> Line <ref>). The collated 𝐱_u and r are provided as an input for PSP (ϕ) (Algo. <ref> Line <ref>). The output of ϕ forms candidate solutions for the LL PS or initial solutions for LL search (Algo. <ref> Line <ref>).




The above generated solutions are evaluated using the LL objective/constraint functions. Then, for the generations where PS approximation is directly used, the ND of the above solutions are returned back to the UL. For the generations where LL search is invoked (i.e., every γ generations; Algo. <ref>, Line <ref>), the ND solutions are used to seed the first LL population (Algo. <ref>, Line <ref>). 








For the LL searches in the first UL generation, the initial population is randomly generated, while subsequently it comes partially from predicted LL ND solutions. If ND solution size is smaller than pre-defined population size n^l,  then the rest of the solutions are filled in with randomly generated solutions (Algo. <ref> Line <ref>). The intent behind this is that if the initial solutions generated by PSP model ϕ are clustered in a small region, these randomly introduced candidate solutions can help maintain diversity to a certain extent. Inferred from vector v, an LL solution can be divided into (𝐱_li^U, 𝐱_li^L), where 𝐱_li^U refers to LL variables that only appear in UL objective. If there is no VAA in the problem (which is most often the case), then 𝐱_li^U is empty. LL search discussed above is only applied to 𝐱_li^L (Algo. <ref> Line. <ref>). This brings benefits of smaller search space and potentially quicker convergence speed to LL search, for the cases where 𝐱_li^U is non-empty.

After first population P_1^l is formed and evaluated, evolutionary operators are used to generate child population C_g (Algo. <ref>, Line <ref>). Environmental selection (Algo. <ref>) is then invoked to select the surviving population. This iterates until stopping condition is met. Then ND solutions from the last population are selected (Algo. <ref>, Line <ref>).
With 𝐱_l^L variables determined in L, we still need to identify the values for 𝐱_l^U by calling Algo <ref>, if VAA is present (i.e., 𝐱_li^U is non-empty), as discussed next.









 §.§ Additional UL search in case of VAA


This module determines the value of 𝐱_l^U part of the search result L from Algo <ref>. For each solution 𝐱_li = (𝐱_li^U, 𝐱_li^L*) in L, to determine values for the part 𝐱_li^U, MOEA search is conducted to optimize the UL objectives with 𝐱_li^U as variables, while 𝐱_li^L* stay fixed. However, running a standard MOEA for every single 𝐱_li will require exorbitant FE consumption. To expedite this, we utilize solution transfer in this MOEA search. As shown in Algo. <ref> Lines <ref>-<ref>, once there are one or more solutions available in L^*, then when new search to be started, solution transfer is using solutions from L^*. For this, the solution 𝐱_lj^L* that has the smallest Euclidean distance from 𝐱_li^L* (i≠ j) is identified. Then, 𝐱_lj^U* corresponding to this solution is inserted into the initial population for the search corresponding to this 𝐱_li^L. Once 𝐱_li^U* is determined for each 𝐱_li^L*, the updated L^* is returned. 








 §.§ Environmental Selection


In the evolutionary process of Algo. <ref> (Line <ref>), and Algo. <ref>, (Line <ref>) new population is selected from the sorted combined parent and offspring population. In order to maintain diversity in the population, we adopt distance based subset selection (DSS) <cit.> in this step when ND front size has more unique solutions than predefined population size. As shown in Algo. <ref>, Line <ref>, from sorted population, we first identify the number N of unique ND solutions in {𝐱}_g. If N is smaller than predefined population size, then new population is selected in order of ranking. If N is larger than predefined population size, then DSS selection is applied to {𝐱}_g's ND front solutions. DSS selects solutions iteratively by choosing the one with the maximum distance to already selected solutions. As in Eq. <ref>, 𝐱_i refers to candidate solution, 𝐱_j^s refers to solutions already selected, and k is the number of selected solutions. In each iteration, the solution that has highest value of d_𝐱_i is selected. 

    d_𝐱_i = min{d(𝐱_i, 𝐱_1^s)), d(𝐱_i, 𝐱_2^s), ...d(𝐱_i, 𝐱_j^s), ...d(𝐱_i, 𝐱_k^s)}







 §.§ Termination condition


To terminate the algorithm, we adopt the condition proposed in <cit.>, which  
attempts to detect the stability (stagnation) of the evolving population based on measurements capturing convergence and diversity. Convergence is measured by the following normalized ideal (z^*) difference and nadir (z^nad) difference. 


    Δ z_t-1, t^* = max_i=1^M z_i^*(t-1) - z_i^*(t)/z_i^nad(t) - z_i^*(t)



    Δ z_t-1, t^nad = max_i=1^M z_i^nad(t-1) - z_i^nad(t)/z_i^nad(t) - z_i^*(t)

If Δ z_t-1, t^* and Δ z_t-1, t^nad do not change significantly in the past ω generations by comparing their maximum metric values within a threshold ϵ, then the algorithm is deemed to have sufficiently converged. Here t refers to a given generation during evolution. M refers to the number of objectives.

The diversity is assessed using the following metric (but as noted in <cit.>, is not entirely independent of convergence). 


    ϕ(t) = IGD(P^t(t-1), P^t(t))
 

IGD refers to inverted generational distance <cit.>. Suppose 0≤ t ≤τ, τ is current generation, the normalized i^th objective value of j^th point in generation t with regard to τ generation is computed as


    P_i^τ, j(t) = P_i^j(t) - z_i^*(τ)/z_i^nad(τ) - z_i^*(τ)


A sliding window ω is used to compute the metrics. If over ω generations, the maximum values of three metrics given below are no greater than a threshold ϵ, the algorithm is terminated. 


    max(Δ z_τ, τ-1^*, ..., Δ z_τ-ω, τ -ω -1^*) ≤ϵ
         max(Δ z_τ, τ-1^nad, ..., Δ z_τ-ω, τ -ω -1^nad) ≤ϵ
        max(Δϕ_τ, τ-1, ..., Δϕ_τ-ω, τ -ω -1) ≤ϵ



Similar to IGD, it is also possible to monitor HV to formulate a termination condition, as done in some of the other works <cit.>. Given that the true PF is not known a priori while solving the problem, either of these can only measure whether the population has stabilized, without guaranteeing convergence to the true optimum. We choose IGD based stopping condition as the default in our algorithm, merely to be consistent with the use of IGD metric later for performance assessment. IGD is a more reliable metric for measuring the final performance for BLMOPs, especially where deceptive problems are involved. For deceptive problems, a higher HV value may not reflect true performance of an algorithm, since it may be an artifact of sub-optimal solutions at the LL. The UL solutions in such cases may appear to dominate the true PF and generate high HV value, contradicting its true performance. IGD metric on the other hand becomes worse when PF approximation is away from the true PF, irrespective of which side it is located, providing a more accurate measurement. We adopt above IGD based stopping condition as default in our studies. However, some exceptions lie in  comparison with state-of-the-art algorithms where we have used termination based on HV or/and number of evaluations to be consistent with that used in other works.




§ EXPERIMENTS AND DISCUSSION



Ten problems are used in our experiments from study <cit.>, which are listed in Table <ref>. Five of them (i.e. TP1, TP2, DS1D, DS2D, DS3D) exhibit some level of deceptiveness. Theoretical PF can be derived for all the test problems to analyze the performance of the algorithms. To generate n points approximately uniformly distributed on the PF, we first over-sample by generating 2n approximately uniform points on PF, utilizing the equations relevant to the PF of the problem. Then we use distance based subset selection <cit.> to select n solutions (n=1025 to ensure near uniformity) as the final PF. For calculating HV when comparing to state-of-the-art algorithms, reference point is set 1.1 times of the maximum objective values of the true PF, without normalization, to be consistent with their settings.




For the evolutionary search part of the proposed framework, differential evolution (DE) <cit.> and polynomial mutation <cit.> are used for generating offspring population. The crossover rate is set to 1 and scaling factor is set to 0.5 for DE operator. As for PM, its mutation probability is set to 1/D (D is the number of variables) and mutation index to 20. The first generation of PSP-BLEMO relies on MOEA for its LL search. This MOEA has fixed generation size set relative higher to 300 to ensure first training data has good quality. Stopping conditions as discussed previously are activated during the remainder of the search. 21 independent runs are conducted for each setting. Wilcoxon rank-sum significant test is used to draw conclusions regarding statistical significance. 

For the NN structure used in PSP-BLEMO, one hidden layer is used. The size of hidden layer nodes is set to two times of the size of input nodes or output nodes, whichever is higher. Two parameters of PSP-BLEMO, re-train gap γ and training data size ds, are determined empirically through experiments, discussed shortly in more detail; however it is to be noted that with sufficient training datasize, algorithm performance is observed to be not too sensitive to γ. Matlab  tool is used to train the NN with all parameters at their default values. For example, learning rate is 0.1, loss function is mean square error between network output  and target output. Levenberg-Marquardt is the default training algorithm. NN training stops if validation error fails to decrease for 6 (default) iterations. The split ratio for training, validation and test is 70%, 15% and 15%, respectively.

In the following subsections, we present three sets of experiments. In the first set, we empirically determine two key parameters of PSP-BLEMO, namely the training gap γ and the data size threshold ds. After their values are fixed, we continue on to the second set of experiments, which compare the performance of PSP-BLEMO with six state-of-the-art algorithms, namely cG-BLEMO <cit.>, MOBEA-DPL <cit.>, SMS-MOBO <cit.>, BLMOCC <cit.>, mf-BLEAQ <cit.>, H-BLEMO <cit.>. For these comparison studies, both UL and LL population sizes are fixed to 20, same as in the most recent study <cit.>. The parameter settings for SMS-MOBO, are consistent with <cit.> too. As discussed in above stopping condition section, for these experiments, HV stopping condition <cit.> is used to be consistent with these methods. The thresholds for HV stopping condition are 1e-3 over ω=10 generations, same as compared methods in <cit.>. For the additional UL search, population size is set to 5 and first search generation size is set to 80. During the UL search assisted by transferred solutions (Algo. <ref>), HV stopping condition is used for consistency.  In the last set of experiments, we compare PSP-BLEMO with its variant where the training is done only once (referred to as one-shot (OS) method) in order to save on evaluations. It reveals some interesting observations are when larger population sizes (proportional to number of variables) are used, unlike the previous two experiments.




 §.§ Empirical experiments on key parameters


In this first experiment, we use empirical analysis to gauge the effect of γ and ds. For γ, we investigate 5 different values, 5, 10, 15, 20 and infinite (Inf). Inf here means that once the data reaches ds size, the model is trained and thereafter never updated. As for ds, we tested 5 values: 5e2, 1e3, 2e3 and 5e3. Test bed involves all 10 problems shown in Table <ref>. For each combination of γ and ds, 21 runs are conducted for each problem. Both UL and LL stops using default IGD stopping condition, with thresholds (ϵ) being 1e-2 over 5 generation window (ω). The statistics on IGD median values and FE consumption are reported in the supplementary (Section 1). Here, we include a summary of statistical significance tests in Table <ref>. For each γ setting, we run Wilcoxon rank sum tests (0.05 significance level) between each data size setting to the others on each problem, and count how many times this setting shows significant better performance. For example, when γ=5 is fixed, run statistical tests on test problems between (1) ds=5e2 vs ds=1e3, (2) ds=5e2 vs ds=2e3, (3) ds=5e2 vs ds=5e3, resulting in 30 tests (10 test problems × 3 pairwise comparisons). Setting ds=5e2 shows significant better, equivalent and worse performance in 0, 13 and 17 tests, respectively. 






The left half of Table <ref> basically evaluates which data size ds setting has most win cases when γ setting is unchanged. It can be seen that ds=5e3 has highest win cases in all γ settings. Therefore, in following experiments, we fix training data size ds=5e3. Following the same approach, the right half of Table <ref> shows the number of win cases for each γ setting when ds is fixed. We can see that when ds=5e3, γ=10 has most win cases. Consequently ds and γ are set to these two values for further experiments. 





 §.§ Comparison with state-of-the-art approaches



The compared algorithms include the most recently proposed cG-BLEMO <cit.> along with MOBEA-DPL <cit.>, SMS-MOBO <cit.>, BLMOCC <cit.>, mf-BLEAQ <cit.> and H-BLEMO <cit.>. For BLMOCC and H-BLEMO we inherit the reported results in their corresponding papers, due to unavailability of the codes. The key idea of PSP-BLEMO is to replace LL PS with predicted value, cG-BLEMO and mf-BLEAQ are the most closely related algorithms for comparison, as both involve modules predicting LL solutions. For mf-BLEAQ, the source code was not available and hence we implemented its LL PS predictor, i.e. set valued mapping (SVM) and embedded it into PSP-BLEMO framework for comparisons.






In Table <ref>, the mean and standard deviation of HV and IGD values obtained by the six algorithms are shown. 
The symbols ↑, ↓ and ≈ show significantly better, worse or equivalent result, respectively, compared to PSP-BLEMO. For BLMOCC and H-BLEMO the statistical tests are conducted through confidence interval test <cit.>. If two algorithms' 95% confident intervals are overlapping, then there is no statistical difference between two algorithms, otherwise, statistical difference is reported. For others, Wilcoxon ranksum test is used. Since our experiments include deceptive problems, HV performance can be misleading (as discussed earlier in Section <ref>). We therefore mainly focus on IGD performance, but still show HV results for completeness. With regards to cG-BLEMO, we observed certain trade-offs between function evaluations and performance, so we have reported two values. The second column in Table <ref> shows the performance of cG-BLEMO after its run is completed. The third column is also cG-BLEMO, but marked with FE truncated (tr). The metric values of cG-BLEMO are extracted when its total FE consumption exceeds total FE consumption of PSP-BLEMO for each run (across 21 runs). The process is as follows. Both PSP-BLEMO and cG-BLEMO runs are sorted based on the FE consumed and then the metric values of cG-BLEMO are extracted when its total FE consumption exceeds the FE consumption of PSP-BLEMO for the corresponding run in the sorted order. This truncation allows us to see performance difference between PSP-BLEMO and cG-BLEMO for similar computational budgets.

In terms of mean IGD (Table <ref>), PSP-BLEMO shows the lowest value for 5 out of 10 problems, while cG-BLEMO has the lowest for 1. In terms of significance test, cG-BLEMO performs better than PSP-BLEMO in 4 problems, and worse or equal in the other 6 problems. PSP-BLEMO shows better mean and standard deviation values especially in deceptive problems when compared with cG-BLEMO. This implies that performance variation of cG-BLEMO is much higher when dealing with deceptive problems. This can be explained as follows. cG-BLEMO inserts UL evaluations on LL solutions during LL search. If its LL prediction module initializes LL solutions that are not close enough to PS, UL objectives of such LL solutions may dominate the true UL PF. 
Consequently, UL search may be misguided by evolving around corresponding UL solutions which show inaccurate UL objective values; resulting in poor search performance. In this respect, mean and standard deviation reflect how well LL prediction maintains its performance. PSP-BLEMO shows relatively stable performance compared to cG-BLEMO. 











Next, we look into FE consumption of all algorithms, reported in Table <ref>. For brevity, raw values are reported for PSP-BLEMO, while for the rest of the algorithms, ratio of their FE to that of PSP-BLEMO is reported. It can be seen that cG-BLEMO uses around 7 to 40 times more FE than PSP-BLEMO at the UL. Likewise, on the LL, cG-BLEMO uses around 3 to 15 times more FEs than PSP-BLEMO. In previous IGD metric based performance, cG-BLEMO showed better statistical performance in 4 out of 10 cases, but it also consumes much higher resource. To have one of these two metrics (IGD or FE) on the same level, we report cG-BLEMO IGD performance when its FE consumption exceeds PSP-BLEMO's in column 3 Table <ref>.  Statistical analysis shows that PSP-BLEMO performs significantly better than cG-BLEMO in 8 out of 10 cases, suggesting the improved performance of PSP-BLEMO for lower evaluation budgets. To note that  DS4 and DS5 are two special cases, where additional UL search as suggested in Algo. <ref> are invoked. It can be observed that PSP-BLEMO consumes much more FEs for such problems especially on the UL. However, even when considering similar FEs, its performance still remain competitive to cG-BLEMO. In supplementary (Section II), we also present a comparison of PSP-BLEMO with cG-BLEMO for deceptive versions of DS4 and DS5 (formulated in <cit.>) to highlight the advantage of PSP-BLEMO for such problems. Additionally, in the supplementary (Section III), we report FE comparisons based on when these algorithms achieve the same level of IGD as BLMOCC <cit.>. This method for comparison was adopted in <cit.>, hence we report them for completeness. PSP-BLEMO is observed to be competitive in these comparisons, consistent with the results above.

The above performance comparison between PSP-BLEMO and cG-BLEMO is also visually reflected in convergence plots of median IGD runs shown in Fig <ref>. Typically, PSP-BLEMO uses much fewer FEs to converge, at a faster rate than cG-BLEMO. Besides FE and IGD, there is also additional information implicit in the plots. For TP1 and DS1D, IGD values of cG-BLEMO decreases first and then rises. This is reflection of the algorithm performance affected negatively by the deceptive problems. In contrast, such case do not occur for PSP-BLEMO. For problems DS4 and DS5 (the cases with VAA), it can be observed that there is suddenly drop in IGD value in the first generation for PSP-BLEMO. It shows the effect of UL additional search to improve the IGD. In Table <ref>, we compare the average runtime of two algorithms. Consistent with FE consumption, it can be seen that PSP-BLEMO uses much less time to converge in 8 out of 10 problems. For DS4 and DS5, the two algorithms use similar time.








SVM-BLEMO also uses LL PS prediction. However, its performance is significantly worse than PSP-BLEMO, as evident from Table <ref>. The IGD values are much higher than rest of the compared methods. This suggests that the LL search results of SVM-BLEMO are relatively further away from the optimal. Due to lower accuracy of LL PS prediction, it is not ideal for SVM based predictor to skip LL search entirely (for up to γ generations) in the proposed framework.

Four other state-of-the-art algorithms, namely BLMOCC, H-BLEMO, MOBEA-DPL and SMS-MOBO, are less similar conceptually to PSP-BLEMO. BLMOCC adopts co-evolution on both levels, while H-BLEMO involves hybridization via local search. One can observe from Table <ref> that PSP-BLEMO has lower mean IGD values than BLMOCC and H-BLEMO. In terms of FE, compared to BLMOCC, PSP-BLEMO uses fewer FE on the LL, and and more FE on the UL. However, the FE consumption on LL is usually 10 to 100 more times than UL (seen from Table <ref>). Thus, overall, PSP-BLEMO holds advantage over BLMOCC in both IGD metric and FE consumption. Similar observations are made with respect to H-BLEMO; additionally, PSP-BLEMO shows fewer evaluations also on the UL. MOBEA-DPL adopts two populations on the LL search, of which one population is evaluated on the UL. Therefore deceptive problems are challenging for MOBEA-DPL.  In terms of FE consumption, it uses up to 30 times more FE than PSP-BLEMO on the UL and up to 14 times more FE on the LL. SMS-MOBO shows competitive performance in 4 problems, but its performance needs up to 100 times more FEs (on the LL) to achieve this performance. SMS-MOBO also faces challenges on deceptive problems, where it uses more FEs but returns rather high IGD values. Compared to SMS-MOBO, PSP-BLEMO shows stable performance across different problem types. 


Lastly, we explicitly observed the quality of the LL PS predictions, to account for the improvements in performance achieved by PSP-BLEMO. In Fig. <ref>, we compared LL prediction results on PSP-BLEMO and cG-BLEMO. In Fig. <ref>, it shows the prediction results obtained by two algorithms after first UL generation. PSP-BLEMO's prediction is not close to PF for this case. The reason is that after the first generation training data accumulated was only 400 (70% training data therefore uses 280 data points). The size of training data is not enough to capture the UL solution to LL PS mapping. cG-BLEMO's PS prediction is relatively better in terms of being closer to PF. However, when the training data size increases to 1000 (70% training is 700 data points), PSP-BLEMO's prediction is much closer to the true LL PF and well distributed as shown in Fig. <ref>. As for cG-BLEMO, its prediction improves, but there is still significant amount of search effort required to move towards LL PF.


To summarize, PSP-BLEMO shows competitive performance in terms of UL IGD, FE consumption, and average runtime when compared with other state-of-the-art algorithms. Especially, when the target problem has deceptive nature and the number of allowable evaluations is relatively low. The competence of  PSP-BLEMO is supported by the LL PS prediction quality, which helps it skip  LL search for most of the generations with low impact on the quality of solutions obtained.




 §.§ Comparison study with `one-shot' PSP-BLEMO


Above comparison with state-of-the-art algorithms shows competitive performance of PSP-BLEMO in terms of accuracy and FE consumption. The experimental settings used in previous studies (e.g. population size 20) have been used for fair and consistent comparisons with the state-of-the-art algorithms. However, in this section, we would like to look into a variant of PSP-BLEMO and its performance when experimental settings and problem settings are varied. PSP-BLEMO, in its default form, updates NN every γ generations. To attempt further savings in FE, one may consider a version of PSP-BLEMO that forgoes re-training entirely and train the model only once, in the first generation. We refer this version one-shot (OS) PSP-BLEMO. During the run, the LL PS is always predicted, never obtained through search. At the end of the run, for the final reporting, the solutions of the last population are evaluated by running predicted LL assisted search on the ND solutions in UL archive. This process helps improve the solution quality on the LL, to more accurately reflect the quality of returned UL solutions. For a fair comparison, such a re-evaluation scheme is also applied to PSP-BLEMO in these experiments. 
In what follows, we use the shorthand OS and PSP to distinguish these two variants. Further, as a baseline, we use nested evolutionary algorithm (NE), where both UL and LL search are carried out by MOEA described earlier in Section <ref> without the use of any predictions. For test problems, we expand the test bed to include problems of different variable sizes, so that the scalability of algorithms can also be observed. The detailed problem variable settings can be found in Table <ref>. 













In the first part of this experiment, we keep the population size fixed to 20, as done previously. The performance in terms of normalized UL IGD, normalized LL IGD and total FE are summarized in Table <ref>. The calculation of normalized LL IGD is as follows. Firstly, for every UL search result, we compute the normalized LL IGD for each UL solution. Next, the mean of these normalized LL IGD values is associated with each search result. Subsequently, we identify the maximum and minimum normalized IGD values across all three algorithms (considering all seeds) to establish the range for normalizing each LL IGD result between 0 and 1. The reported normalized LL IGD is the final result from above 3 steps.
Statistical test results are shown in symbols ↑_i, ↓_i and ≈_i. This i means compared to (i) algorithm, current algorithm (on the left of symbols) performs better, worse or equivalent. For FE consumption, except for DS4 and DS5, it is not surprising that OS is the best, and PSP uses more FE than OS, but both of them use fewer FE than the baseline.
The larger the number of variables, the more savings on FE the proposed method can achieve. 
For DS4 and DS5, due to additional search on the UL, PSP has more FE consumed than baseline, but this number decreases as the number of variables increases as shown in the last section of Table <ref>. For OS, its FE is higher than baseline when variable size is small (S1), but in the rest of the settings, its FE consumption becomes much smaller than baseline and PSP.


It can be noted that PSP shows better performance in all variable settings. For S1, PSP performs equal to or better than baseline in all problems on the UL, 9 out of 10 on the LL. For S2 and S3, PSP performs equal or better than baseline  in all problems on UL and LL. Between PSP and OS, similar trend can be observed on the UL. 
For all problems except one across three settings, PSP performs better than or equal to OS. 
LL situation varies mainly in S3 section, where where OS outperforms PSP 4 out of 10 cases on LL. But when we check the corresponding UL performance, e.g. DS2, DS3 and DS2D etc, OS has much worse UL IGD values. Overall, PSP presents most consistent performance in this setting. The performance of OS is also inline with expectation, i.e. with a population size of 20 and for small number of variables, it would have accumulated enough training data for appropriate prediction accuracy. However, for larger number of variables, its prediction accuracy is compromised affecting the overall performance.


In the second part of this experiment, we set the population size to be 10D and maximum generation size to be 10D (while also retaining the IGD termination criterion), making them proportional to the problem variable size. D refers to corresponding level variable size. For UL, it means D_u, while for LL it means D_l. Table <ref> shows results of these experiments. Again, it is seen that OS and PSP use much fewer FE compared to the baseline except DS4 and DS5. For DS4 and DS5, similar to the first experiment, when problem variable size increases, this extra FE consumption also decreases compared to those used by baseline. Between PSP and baseline, on both levels and across all problem settings, PSP performs significantly better.  
Between PSP and OS, the performance is not as overwhelmingly in favor of PSP as it has been in the previous experiments.


For S1 and S2, PSP still shows better performance on both levels, while for S3, PSP performs better than OS at UL in 4 problems and worse in 3 problems. On the LL, PSP performs better or equal to OS in all cases. For S3, the population size and generations are significantly larger in this case (e.g. 100 for 10 variables). Therefore, OS has same number of training points as PSP, except that PSP can also update the models later. The competitive performance of OS suggests that the proposed LL PS prediction model, being built on a larger dataset, can maintain its performance quite well over generations. This observation actually also resonates with the Table <ref> last column about test results on γ when ds=5e3. Although according to wining cases, we choose γ=10 as our parameter value, for majority of the cases, the differences are not significant (i.e. majority of the cases are `equivalent'), γ=10 is only marginally better than others. So, the OS performance can be competitive when ds is large. Overall, we still tend to choose PSP with re-training (every γ generations) as the default. Although the FE consumption on DS4 and DS5 is high, it is still competitive with the state-of-the-art, and such problems are not too common in BLMOPs. Considering the performance at both UL and LL, it is the most consistent and competitive strategy across different settings and provides an opportunity for model correction intermittently. But the insights gained from this last set of experiments provides an avenue for further FE reduction. The use of OS strategy can be beneficial towards this endeavor without significant compromise on performance, provided a sufficiently large dataset is used to train the models.




§ CONCLUSION AND FUTURE WORK


Solving BLMOPs using a nested approach requires exorbitant number of function evaluations. Therefore, it is of interest to develop techniques that can reduce the computational expense while achieving good solution quality. In this paper, we develop a Pareto set prediction assisted search (PSP-BLEMO) towards addressing this challenge. To achieve this, the archive of evaluated solutions is transformed using a helper variable, and an NN model is trained on the resulting dataset. The model is then used to predict the PS approximation, used to seed the LL population or bypass the LL search entirely. The model itself can be re-trained periodically using additional UL solutions that undergo LL search every γ generations. For saving additional FEs, a `one-shot' approach can be applied where the model is not updated during the run.

Extensive numerical experiments were conducted to test the efficacy of the proposed approach in terms of, solution quality, FE consumption and runtime. Key parameters and variants of the algorithm were also analyzed. Comparisons with the state-of-the-art algorithms shows competitive and consistent performance of PSP-BLEMO. The performance is particularly notable where deceptive functions are involved. Moreover, it was concluded that the performance of the algorithm was not too sensitive to the training gap γ, as long as there is sufficient training data available to build accurate PSP models.

In the future work, further improvements in the algorithmic performance will be investigated, e.g., by enhancing the modeling through use of more carefully selected subset from the data. Currently, the most recent LL search results are used as training data. Its distribution might not be uniformly distributed over the entire PF. Decomposition based segmentation on objective space can provide opportunity in improving this data quality. 
Moreover, extensions for higher number of objectives could also be considered. The way in which helper variable r is used in this study is particularly suited for two-objective problems, since arranging r and f_1 corresponds to sequentially traversing the span of the PF. However, for problems with more objectives, the composition of r will need to be further refined. One potential way is to utilize reference vectors obtained through systematically sampled points using normalized boundary intersection method, commonly utilized in decomposition-based MOEAs. In the current state of the field, more than two objectives are rare, hence this would be an interesting extension of this study in the future. Lastly, extension of some other recent modeling structures, such as hypernetworks, to solve BLMOPs could also be explored.



ACM-Reference-Format




