Solving combinatorial bi-level optimization problems using multiple populations and migration schemes

Abstract

In many decision making cases, we may have a hierarchical situation between different optimization tasks. For instance, in production scheduling, the evaluation of the tasks assignment to a machine requires the determination of their optimal sequencing on this machine. Such situation is usually modeled as a Bi-Level Optimization Problem (BLOP). The latter consists in optimizing an upper-level (a leader) task, while having a lower-level (a follower) optimization task as a constraint. In this way, the evaluation of any upper-level solution requires finding its corresponding lower-level (near) optimal solution, which makes BLOP resolution very computationally costly. Evolutionary Algorithms (EAs) have proven their strength in solving BLOPs due to their insensitivity to the mathematical features of the objective functions such as non-linearity, non-differentiability, and high dimensionality. Moreover, EAs that are based on approximation techniques have proven their strength in solving BLOPs. Nevertheless, their application has been restricted to the continuous case as most approaches are based on approximating the lower-level optimum using classical mathematical programming and machine learning techniques. Motivated by this observation, we tackle in this paper the discrete case by proposing a Co-Evolutionary Migration-Based Algorithm, called CEMBA, that uses two populations in each level and a migration scheme; with the aim to considerably minimize the number of Function Evaluations (FEs) while ensuring good convergence towards the global optimum of the upper-level. CEMBA has been validated on a set of bi-level combinatorial production-distribution planning benchmark instances. The statistical analysis of the obtained results shows the effectiveness and efficiency of CEMBA when compared to existing state-of-the-art combinatorial bi-level EAs.

Appendices

Appendix 1: A comparison between the Das-&-Dennis and the DSDM methods.

This first appendix illustrates the difference between the Das-&-Dennis (Das and Dennis 1998), and the DSDM (Chaabani et al. 2015) methods as shown in Fig. 6. For the case of a continuous search space, the Das-&-Dennis method could be used. However, for the case of a discrete search space, the Das-&-Dennis method is inapplicable. The DSDM method is a variant of Das-&-Dennis method and it could be used to generate a set of points in discrete search spaces. We must mention here that the Das-&-Dennis method generates a set of solutions in the objective space, while the DSDM method works in the decision space. The distribution of the reference points for the Das-&-Dennis method with 3-objective optimization problem (M = 3) and a spacing of $\delta = 0.2$ (P = 5) is presented in Fig. 6a. Thus, 21 reference points (H = 21) are generated in a normalized hyper-plane. We mention here that the reference directions are represented by the lines, which are constructed from the origin to each of these reference points. Figure 6b illustrates the obtained results for the DSDM method with three decision variables where the domains are: $D_{x_{1}}$= [0,2,5,13], $D_{x_{2}}$= [4,7,9,17], and $D_{x_{3}}$= [5,8,11,16]. In order to generate the reference points with the DSDM method in a discrete space, a uniform spacing noted $\delta _{i}$ is calculated for each decision variable as follow: $\delta _{i} = max_{i}/P$ (i is the decision variable number, P is a fixed parameter based on the dimension of the problem). Thus, the obtained $\delta _{i}$ for a P = 3 are: $\delta _{1} = 13/3=4$, $\delta _{2} = 17/3=5$, and $\delta _{3} = 16/5=5$. After that, the range values ($R_{i}$) are generated for each decision variable by adding the first value in each range on a set $R_{i}$, and determining the following $R_{i}$ members that obey the $\delta _{i}$. Thus, the obtained range values are as follows: $R_{1} = [0,5,13]$, $R_{2} = [4,9,17]$, and $R_{3} = [5,11,16]$. The obtained solutions for this example are: (0,4,5), (0,9,11), (0,17,16), (5,4,5), (5,9,11), and (13,4,5).

Fig. 6

Illustration of the difference between Das-&-Dennis method, and DSDM method

Appendix 2: Illustration of the interaction process between the two levels in a BLOP

This second appendix is devoted to illustrate in details the interaction process between the upper-level and the lower one in bi-level optimization. Indeed, these two levels are dependent on each other because: (1) the upper-level solution vector $x = (x_{u}, x_{l})$ could not be evaluated without optimizing its corresponding solution subvector $x_{l}$ and (2) the lower-level solution $x_{l}$ could not be evaluated without receiving its corresponding upper-level solution subvector $x_{u}$ as fixed parameter (Huang and Wang 2020; Sinha et al. 2020, 2017b). For this reason, to precisely evaluate any upper-level solution $x = (x_{u}, x_{l})$, we need to effectively approximate the optimal lower-level solution $x_{l}^{*}$ that corresponds to the upper-level solution subvector $x_{u}$, which acts as a fixed parameter within the lower-level search process. In the following, we describe the sequence of steps required to well compute the fitness function of an upper-level solution $x = (x_{u}, x_{l})$ as illustrated by Fig. 7. First, the subvector of upper-level decision variables $x_{u}$ is passed as a fixed parameter to the lower-level search algorithm. Second, this latter evolves a population of lower-level solutions (i.e., lower-level decision variable vectors) for a number of generations with the aim to approximate the optimal follower solution $x_{l}^{*}$ corresponding to the parameter $x_{u}$. Third, once the lower-level algorithm termination criterion is met, the obtained approximation of $x_{l}^{*}$ is sent to the upper-level and thus the upper-level solution vector becomes $x = (x_{u}, x_{l}^{*})$. As illustrated by Fig. 7, to evaluate the upper-level solution $x^{A} = (x_{u}^{A}, x_{l}^{A})$, $x_{u}^{A}$ is passed as fixed parameter to the lower-level EA. The latter executes an evolutionary process to approximate the optimal lower-level solution ${x_{l}^{A}}^{*}$ corresponding to $x_{u}^{A}$. Based, on Fig. 7, the obtained approximation is the solution vector E. Finally, E is assigned to $x_{l}^{A}$; and thus the upper-level solution becomes $x^{A} = (x_{u}^{A}, E)$. The fitness function of $x^{A}$ could now be computed. It is important to note that the better the approximation E is, the higher the precision of the fitness evaluation of $x^{A} = (x_{u}^{A}, x_{l}^{A})$ is. We conclude that the quality of any upper-level solution vector $x = (x_{u}, x_{l})$ depends on two main factors that are: (1) the quality of the $x_{u}$ sub-vector that depends on the values of its components (upper-level variables) and (2) the quality of its obtained corresponding optimal lower-level $x_{l}$ solution approximation (lower-level variables). For this reason, the lower-level algorithm plays a crucial role in the precise evaluation of each upper-level solution quality (upper fitness value). Differently speaking, the better the quality of $x_{l}^{*}$ for a particular upper-level solution $x = (x_{u}, x_{l})$ is, the higher the precision of the fitness computation of $x = (x_{u}, x_{l})$. This high precision allows the bi-level algorithm to compute the fitness function values of the upper population with more exactitude and thus detecting promising search directions in the upper-level search space more effectively and efficiently. This is the key characteristic of the interaction process between the upper-level algorithm and the lower-level one in bi-level optimization. It is important to note that if we consider the case of a poor lower-level search process, the upper level fitness values computations will be biased with noise; which will significantly deteriorate the search process of the upper-level algorithm. For an extremely bad case, a very poor lower level search process could make the upper-level algorithm behavior acting like a random search. For this reason, there is a need to evaluate the performance of the compared algorithms’ lower-levels search methods. To do so, Legillon et al. (2012) proposed the direct rationality and the weighted one (defined in Sect. 4.3) as two metrics to evaluate the performance of a lower-level algorithm (cf. Algorithm 4). The former expresses the mean number of times the lower-level algorithm was able to improve the solutions with respect to the population of the previous generation, while the latter evaluates the mean quantity (magnitude) of the improvements. In summary, we conclude that: (1) the lower-level optimality approximation is considered in CEMBA as a necessary constraint that appears as one of the upper-level constraints and (2) the lower-level performance should be evaluated for all compared algorithms to assess the quality of its contribution in the precise evaluation of the upper-level solutions, and hence in the effective exploration of promising search directions in the upper-level search space.

Fig. 7

A general sketch of two upper-level solutions evaluations (A and B) in a BLOP

To empirically show the importance of the optimization of the lower-level objective function (while respecting its constraints), we have conducted a set of experiments using two versions of each of the two following algorithms: CEMBA and CODBA-CRO:

• 1. CEMBA-WLO and CODBA-CRO-WLO are two variants of CEMBA and CODBA-CRO (respectively) that do not optimize the lower-level objective function (WLO means Without Lower-level Optimization); and

• 2. The original CEMBA and CODBA-CRO, which already optimize the lower-level objective function.

The obtained results on the 23 bip instances and the 10 bipr ones are presented in Table 8. We observe that the results of the original CEMBA and CODBA-CRO are extremely far better than those of CEMBA-WLO and CODBA-CRO-WLO, for all test instances. This could be explained by the importance of the optimization of the lower-level objective function. Indeed, when this latter is not optimized (only the upper-level objective function is optimized while respecting all inequality and equality constraints of both levels), the upper-level fitness function computation will not be precise at all. More specifically, based on the illustrative example of Fig. 7, the subvector $x_{l}^{A}$ of the upper-level solution $x^{A}$ will be assigned a lower-level solution that is much poorer than E; which means that it is far from the corresponding follower optimum. In this way, the solution $x^{A}$ will not respect the constraint of the optimality of the lower-level objective function, which will induce an imprecise fitness computation for the upper-level solution $x^{A} = (x_{u}^{A}, x_{l}^{A})$. By doing so for all upper-level population members, the EA will no more be able to guide the search towards the leader global optimum and its behavior will be equivalent to a random search in the upper-level search space, which is the case of CEMBA-WLO and CODBA-CRO-WLO.