Introduction

Multi-objective discrete optimization problems (MODOPs) [1] are the operation optimization problems with discrete variables and multiple objectives, which together with the multi-objective continuous optimization problems constitute the field of multi-objective optimization problems (MOPs) [2]. Compared with continuous MOPs, there is relatively little research on MODOPs although these problems exist widely in real-world applications, such as 0–1 knapsack problems [3], production planning [4], vehicle routing [5], location selection [6], project scheduling [7], resource allocation [8]. Since most of these engineering application problems are subject to constraints, this can be categorized into constrained multi-objective discrete optimization problems (CMODOPs). As a type of optimization problem, it is obvious that the key to solving such problems lies in searching algorithms and constraint processing methods.

Similar to the classic classification of optimization algorithms for continuous MOPs [9], the searching algorithms of CMODOPs can also be divided into the following four categories according to the principle of selecting solutions: dominance-based [10], decomposition-based [11], indicator-based [12], and hybrid selection method-based [13]. Among them, this kind of dominance-based algorithms including NSGA-II and SPEA2 are the most commonly used in CMODOPs [14, 15]. Nevertheless, due to the randomness of the evolutionary algorithm (EA), no matter which kind of algorithm requires a mass of evaluation to obtain satisfying candidate solutions. The evaluation of some practical engineering problems such as antenna design [16], blast optimization [17], trauma system design [18], and power system design [19] are time-consuming and costly so that the expense of obtaining the optimal solution is unaffordable. Therefore, as an efficient tool for expensive optimization problems, the surrogate model [20, 21] has attracted much attention from researchers in different fields. In the expensive optimization, the evaluation of objectives and constraints can only be performed through surrogates trained by data collected from physical experiments, numerical simulations, or daily life. Such optimization algorithms are named as data-driven optimization algorithms (DDOAs) [22, 23], and it is also commonly known as surrogate-assisted evolutionary algorithms (SAEAs). Many existing regression models have been regarded as surrogate models, such as support vector machines [24], radial basis function [25,26,27,28,29], polynomial regression [30], and Kriging model [31, 32]. As far as expensive CMODOPs are concerned, there are two main challenges. One is how to choose and manage the surrogate models for discrete variables, so that the higher approximate accuracy can be obtained at a reasonable calculation cost. The other is how to strike a balance between the optimal objective values and the feasible solutions during the optimization, thus ensuring the quality of the candidate solutions.

At present, most of SAEAs are presented to solve expensive continuous optimization problems, and few scholars also developed such models for expensive CMODOPs. For example, to mitigate the computational expense of models with both continuous and discrete parameters, Swiler et al. [33] investigated and analyzed the performance of three types of response surfaces developed for mixed-variable models. Nguyen et al. [34] proposed a new surrogate assisted genetic programming in terms of automated design problems of dispatching rules for production systems to improve the quality of the evolved rules without significant computational costs. To find the optimal the resource allocation strategy, Wang et al. [35] employed support vector regression surrogate model instead of the time-consuming simulations to show the solution robustness. Aiming at a bi-objective mixed-integer optimization model, Daniel [36] introduced surrogates including radial basis function and polynomial regression to screen the good solutions. The similar problem has been optimized in literature [37] by established the Kriging surrogate model to reduce the identification load of the pollution sources. However, the continuous surrogate model used to solve CMODOPs may generate large approximation errors, which will affect the evaluation precision. Towards this end, Wang et al. [38] came up with a new algorithm called the random forest assisted constrained multi-objective combinatorial optimization (RFCMOCO), in which the tree structure of the random forest (RF) [39] is well suited for approximating problems with discrete decision variables. Nevertheless, training and updating decision trees require considerable computational costs. Therefore, it is necessary to cooperate with the constraint processing strategy with high efficiency, so as to improve the quality of candidate solutions without large computational costs.

Generally speaking, the simplest and most commonly used constraint processing method is the penalty function method [40], which transforms the problems into unconstrained optimization problems by adding penalty factors. However, it is difficult for the user to define parameter values explicitly. To avoid adjusting penalty coefficients, scholars presented a series of new methods named the separation of objectives and constraints, such as stochastic ranking (SR) [41], epsilon constraint handling method [42], and constraint dominance principle (CDP) [43]. Among the successful approaches, the SR method introduced by Runarsson and Yao to strike a balance between better objectives and feasible solutions, which has been widely used in constrained optimization. For example, Balande et al. [44] combined stochastic ranking with an improved firefly algorithm for constrained engineering design optimization problems. Ali et al. [45] integrated stochastic ranking with an evolutionary algorithm to maintain a balance between penalty and objective functions, which can only address constrained continuous problems with a single objective. Besides, to minimize energy consumption, Hector et al. [46] used alternately stochastic ranking, feasibility rules, and 8-constrained methods to analyze the design changes, and validated that stochastic ranking method always obtains better solutions in all proposed torque limits.

Given the simplicity and efficiency of stochastic ranking, it is very suitable to cooperate with a high-accuracy random forest model to optimize CMODOPs, which can effectively improve the calculation speed while ensuring the quality of solutions. In general, fast non-dominated sorting [47] is usually deemed as a sorting basis in the stochastic ranking strategy to handle constrained problems with multiple objectives. However, when faced with such problems with scattered solution space like CMODOPs, the selection pressure of this method will increase, and the diversity of solutions cannot be weakened, which may ultimately result in the poor quality of the selected individuals and waste of computing resources. Considering the fitness assessment from SPEA2 can take into account location information of the population, it may effectively overcome the drawback. Motivated by this consideration, a random forest assisted multi-objective optimization with improved stochastic ranking (RFMOISR) is proposed in this paper for expensive constrained multi-objective discrete optimization problems. More specifically, the main contributions of this work can be highlighted as follows:

  1. (1)

    The improved stochastic ranking strategy is presented by developing fitness mechanism as the novel criterion to ensure the full exploration of solution space and enhance the quality of solutions for expensive CMODOPs;

  2. (2)

    Random forest surrogate models with high precision are incorporated with the improved stochastic ranking strategy to precisely guide the search direction, thereby saving the calculation resources;

  3. (3)

    A flexible probability operator is adaptively adjusted depending on the bias of the different search stages to advance optimization efficiency, which can further improve the algorithm convergence.

The remainder of the paper is organized as follows. “Relative techniques” briefly reviews the relative techniques, while the proposed approach is discussed in “The proposed algorithm”. Numerical experiments are presented in “Experimental studies”, followed by conclusions and future research directions in “Conclusion and future work”.

Relative techniques

Constrained multi-objective discrete optimization problems

Constrained multi-objective discrete optimization problems can be described as:

$$ \begin{gathered} \max F({\mathbf{x}}) = ({{f}}_{{1}} ({\mathbf{x}}), \ldots ,{{f}}_{{\text{m}}} ({\mathbf{x}}))^{{\text{T}}} , \, {\mathbf{x}} \in \,{{D}} \hfill \\ {\text{s.t.}}\;{{ G}}({\mathbf{x}}) = ({{g}}_{{1}} ({\mathbf{x}}), \ldots ,{{g}}_{{\text{c}}} ({\mathbf{x}}))^{{\text{T}}} \le 0. \hfill \\ \end{gathered} $$
(1)

Among them, x = (x1, x2, …, xn)T is a n-dimensional candidate solution, D ⊆ Rn is a n-dimensional bounded discrete decision space. F(x) is an objective vector containing m (m ≥ 2) objective functions, and G(x) is a constraint vector with c constraints. Since multiple objectives conflict with each other, it is difficult to optimize them simultaneously. Therefore, the concept of Pareto dominance was introduced into CMODOPs to obtain an acceptable Pareto front (PF) [48] rather than single-objective optimality. In this work, we focus on solving constrained multi-objective discrete optimization problems with 2–3 objectives.

Random forest

Random forest (RF) [39] is an ensemble machine learning method proposed by Leo Breiman. The main idea of ensemble learning is to build and combine multiple base learners to fulfill better generalization capability. For RF model, multiple sample subsets are first formed by bootstrap sampling. Then, the decision trees are modeled according to the above subsets and the prediction results are output respectively. The mean of the output of all decision trees is generally regarded as the final prediction result. Meanwhile, RF is a simple and supervised method, which is fast and robust to the noise of the target data. Due to the use of bootstrap aggregation or bagging, the random forest model can decrease the error in classification and regression tasks.

The most obvious feature of CMODOPs is that they have discrete decision variables, thus it is a wise choice to employ RF to approximate m objectives and c constraints of such problems. Moreover, a great many theoretical and empirical studies [49] have proved that the random forest algorithm has high prediction accuracy.

Stochastic ranking strategy

Stochastic ranking strategy [44] is the constraint processing approach using bubble sorting technique, which can strike a balance between optimal solutions and feasible solutions. Specifically, two criteria are used to rank the population separately. The first criterion is ranking them as an unconstrained problem with m + c objectives, and the other is ranking them as a constrained problem with m objectives and c constraints. Individuals with better objective values or low constraint violations are selected and retained with a certain probability.

In the SR strategy, the transition from constrained MOPs to unconstrained MOPs with additional objectives can be clearly depicted by Eq. 2, which can significantly improve the performance of constraint processing in SAEAs [50]. Note that such objective and constraint values here are both predicted by surrogates. Furthermore, the fast non-dominated sorting method is often used as a basis for ranking in MOPs.

$$ F({\mathbf{x}}) = \max (f_{1} ({\mathbf{x}}), \ldots ,f_{{\text{m}}} ({\mathbf{x}}),g_{1} ({\mathbf{x}}), \ldots ,g_{{\text{h}}} ({\mathbf{x}}))^{{\text{T}}} . $$
(2)

The proposed algorithm

Overall framework of RFMOISR algorithm

In this section, a novel surrogate-assisted multi-objective evolutionary algorithm is presented to improve the convergence and diversity for expensive CMODOPs under the limited computational budget. The framework of the proposed RFMOISR is outlined in Algorithm 1. More specifically, the algorithm starts building the random forest models with training data for each objective and constraint, and construct the logistic regression model [51] before optimization. It is worth noting that the logistic regression model is used to rectify the boundaries of the feasible regions during optimization to reduce the impact of approximation errors on the constraints as described in the literature [38]. Afterward, the initial population randomly generated is approximated by the surrogates and then the optimization iteration begins. When selecting the parent population for the next generation, an improved stochastic ranking strategy was introduced to strike a balance between satisfying constraints and better objectives. What is more, to advance approximation accuracy, a flexible model management approach is employed to update the surrogate models after individual approximation errors are corrected. The optimization will not stop until the allowed computing resources are exhausted, and then the current optimal solution for CMODOPs will be output.

figure a

Surrogate modeling

The modeling process of the random forest used in this paper is shown in Fig. 1. First, 1000 training samples are randomly generated to form the original training set N, and exact function evaluations are employed to compute the true values of objectives and constraints. Regardless of whether individuals in the training set are feasible solutions, all the true values are put into the RF trainer. Subsequently, the bootstrap method is introduced to resample from N to generate k training subsets randomly. Then, a classification and regression tree (CART) [52] with binary tree structure is constructed in each subset, and according to the minimum square error criterion, the optimal feature is selected for node splitting. Finally, all CARTs are combined into an RF model, and the final output is the average of the k tree outputs.

Fig. 1
figure 1

The generation process of the random forest

Full size image

Each CART branches recursively and stops growing once the segmentation termination condition is reached. In empirical studies, the threshold is set to le−4 * σ2 (σ2 is the variance of all data before the CART is grown). Note that since multi-objective CMODOPs has m objectives and c constraints, the above process needs to be performed m + c times to build the surrogate for each objective or constraint, respectively.

Constraint processing

Considering that the surrogate model introduces some uncertainty, the improved stochastic ranking strategy is presented to deal with the constraints in this part, which not only provide a good guide to the selection of parent population, but also enhance the optimization efficiency. To be specific, the improved stochastic ranking strategy mainly includes two aspects. One is the introduction of the fitness mechanism, which takes into account the information between dominated and non-dominated individuals in light of SPEA2 [53]. Furthermore, a k-nearest neighbor method is added to assign individual fitness comprehensively. As shown in Eq. (3), the fitness values including dominance relations and density information individual i is calculated:

$$ F(i) = R(i) + D(i). $$
(3)

The formulas for each component are as follows:

$$ R(i) = \sum\limits_{{j \in P + Q,{\text{j}} \succ {\text{i}}}} {S({\text{j}})} , $$
(4)
$$ S(i) = \left| {\left\{ {j\left| {j \in P + Q \wedge {\text{i}} \succ {\text{j}}} \right.} \right\}} \right|, $$
(5)
$$ D(i) = \frac{1}{{\sigma_{i}^{k} + 2}}, $$
(6)
$$ {\text{k}} = \sqrt {\left| P \right| + \left| Q \right|} , $$
(7)

where \(\sigma_{i}^{k}\) is the Euclidean distance from individual i to kth neighbor. Hence, the Euclidean distance between the individual i and other individuals from the population P and the external archive set Q is counted and arranged in ascending order. Moreover, the niche method and elite retention mechanism were employed to expand the search range and augment the diversity of the population.

Meanwhile, Fig. 2 visualizes the superiority of the stochastic ranking with the fitness mechanism compared to the fast non-dominated sorting method when the solution space is scattered. Please note that the objectives are to be maximized in Fig. 2. From Fig. 2a there is no obvious distinction between individuals of the same level according to the fast non-dominated sorting method. As a result, the SR strategy may ignore potential individuals, especially those distributed in sparse areas. Whereas as illustrated in Fig. 2b [53], capacity for each individual can be accurately assessed through fitness mechanisms. Besides, a neighbor mechanism is utilized to maintain population distribution. Thus, SR based on this criterion can provide excellent trade-off ability between convergence and diversity.

Fig. 2
figure 2

Comparison of the two different hierarchical schemes for a maximization problem with two objectives f1 and f2

Full size image

The second aspect is the probability operator. In the original SR, a fixed probability is defined to determine whether the comparison with neighboring individuals is based on objectives or constraints. However, EAs have different emphases in different search stages. Therefore, to speed up convergence and save computing costs, we propose a new probability operator of SR. As we all know, the initial population is relatively uniform, and the algorithm is expected to search as wide as possible. With the optimization proceeding, the population gradually converges near the Pareto front, and the algorithm is required to search within the feasible region. Therefore, the fixed probability is replaced to a dynamic probability that can be adaptively adjusted according to search stages, as shown below:

$$ P = {\text{fe}} \div {\text{FE}} \times P_{0} , $$
(8)

where FE is the maximum number of expensive fitness evaluations allowed, fe is the current number of exact evaluations, and P0 is the default probability. The improved strategy is consistent with the search law of EAs, which can further improve the quality of the solutions and the optimization efficiency. Algorithm 2 presents a pseudo code of the improved stochastic ranking (ISR) including Adaptive Selection (AS) and Fitness Mechanism (FM).

figure b

Model managing

In SAEAs, model management is critical to the effectiveness of the surrogate, affecting the optimization performance of the algorithm. To reduce the approximate error, an individual-based model management strategy [54] was employed in this work to correct and update the model in each generation. The approach focuses on the selection of individuals demanding evaluation for each generation. During the optimization, all individuals are evaluated approximately through surrogates. Then, these individuals are sorted according to the fitness mechanism. Based on this order, the potential of individuals can be identified. And the root-mean-square error (RMSEs) between the predicted and true values of each objective for non-dominated individuals are calculated and recorded as E. Next, we use them to modify the predicted objective values F(x) for each individual. Specifically, in the maximization case, the lower bound on the objectives of revised individuals can be expressed as F(x) + E. If the modified individual dominates the non-dominated solutions in the training set, it is regarded as new data to put into the training set for updating the model and ultimately stored into an archive for the next generation. The detail of this strategy is shown in Algorithm 3.

figure c

Experimental studies

To examine the efficiency of RFMOISR, numerical experiments are conducted to make comparisons between RFMOISR and four popular compared algorithms in CMODOPs including RFCMOCO [38], SPEA2 [53], NSGA-II [55], MOEA/D [56].

Testing problems

The CMODOPs are designed to obtain the optimal order or grouping of events or items. Hence, for testing the popularity of the proposed algorithm in solving such problems, common and representative 0–1 multi-objective knapsack problem (MOKP) were selected as benchmark problems. And we assume that the exact evaluation may be computationally expensive. The distinct reason for this is that the computational expense can be easily measured in the problem by limiting the number of exact evaluations. A MOKP problem with m objectives and n items can be expressed as follows:

$$ \begin{gathered} \max f_{j} = \sum\limits_{i = 1}^{n} {v_{i}^{j} x_{i} ,x_{i} \in \{ 0,1\} } ,1 \le j \le m \hfill \\ {\text{s.t.}}\; \, \sum\limits_{i = 1}^{n} {w_{i} x_{i} \le W} , \hfill \\ \end{gathered} $$
(9)

where \(v_{i}^{j}\) is the value of the jth objective corresponding to the ith item, wi is the weight of the i-th item, and W is the maximum load of the backpack. The test questions used in this paper include two parts. One is ten small MOKPs (S-MOKPs) with 10–100 items and 2–3 objectives, which can be obtained by the method of literature [57]. The other is that six widely used large MOKP (L-MOKPs) [58] are considered here, with 250–750 items and 2–3 objectives.

Performance metrics

The inverted generational distance (IGD) [59] is deployed as the index for evaluating the quality of the obtained solution set in terms of both convergence and distribution. The smaller the IGD value, the better the quality of solution set obtained by the algorithm. To compute IGD, a reference set needs to be specified that is a close representation of the Pareto front. For such test problems, we introduce the method from the literature [38] to generate the reference set by 30 runs of NSGA-II (with a population size of 100 runs for 200 generations using exact function evaluations).

The hypervolume (HV) [58] is also adopted as the performance metric in the comparisons to make empirical comparisons between the results obtained by each algorithm. This metric measures the volume of the objective space between the obtained solutions set and a predefined reference point. Since the reference point plays an important role in the calculation of HV, the points far from the objective values are selected as reference points to maintain the accuracy of HV. Given a reference point, a larger HV value means better quality.

The above two metrics evaluate the comprehensive performance of the algorithm, namely convergence and diversity. In principle, the expensive CMODOPs focus more on finding better feasible solutions. Therefore, this paper utilizes the maximum error rate (ME) [60] and generational distance (GD) [61] to measure the convergence of the algorithm. Moreover, a smaller value of ME or GD is preferable, which indicates that the obtained non-dominated solutions have excellent convergence.

Parameter setting

In the experiments, the common parameters, such as the population size (Np) and the maximum exact evaluations (EFmax), were set to the same values in compared algorithms. Furthermore, algorithms use two-point crossover with a probability of 1 and point mutation method with a probability of 0.4. The detail parameter settings in the experiments are given in Table 1. All algorithms are not stopped until the maximum number of exact evaluations is reached, which are performed 30 independent runs on each test problem to obtain a robust optimizing result. All programs were coded in Matlab 2018a on an Intel Core i7 with a 1.99 GHz CPU.

Table 1 Parameter setting of RFMOISR
Full size table

Comparative experiment of ISR strategy

To validate the performance of the proposed constraint handling strategy in the RF model, we compare the ISR strategy with Original Stochastic Ranking (OSR) [41], Stochastic Ranking with Adaptive Selection (SR-AS), and Stochastic Ranking with Fitness Mechanism (SR-FM). The experiment was carried out in the MOKP instances with 30–50 decision variables and 2 or 3 objectives to show the differences between these strategies. And the parent population retained by the NSGA-II using exact function evaluations is considered as the reference set to measure the selection accuracy of the above four strategies. Specifically, the RF models are first trained using 1000 random samples and then used to approximate the objectives and constraints of all the individuals generated by NSGA-II. Subsequently, the individual fitness values are calculated based on two different criteria and regarded as the basis for the individual ascending order. Finally, the SR strategies are employed to select the parent population, and then their selection accuracy is obtained respectively through the reference set. The experiment will be terminated when the maximum number of allowed evaluations is exhausted, of which 1000 exact evaluations are used for building the RF models.

Figures 3 and 4 depict the selection accuracy of these four strategies on such MOKP instances with two and three objectives, separately. It can be seen from Fig. 3a to b that the selection accuracy of ISR can achieve the highest 91% and 92%, respectively, on test instances with 2 objectives, 30 and 50 decision variables, which is significantly better than the other three strategies. Moreover, the performance of SR-FM is second only to ISR, which also shows that a stochastic ranking strategy with an adaptive mechanism can accurately assess the potential of each individual, thereby improving selection efficiency.

Fig. 3
figure 3

The accuracy of selection for different strategies on MOKPs with 30–50 decision variables and 2 objectives

Full size image
Fig. 4
figure 4

The accuracy of selection for different strategies on MOKPs with 30–50 decision variables and 3 objectives

Full size image

What is more, a similar phenomenon can be observed in Fig. 4. From Fig. 4a, the selection accuracy of ISR outperforms the others on the MOKP instance with 30 decision variables and 3 objectives, further demonstrating the superiority of the strategy. In addition, compared with Fig. 4a, the performance of ISR is more prominent on the MOKP instance with 50 decision variables and 3 objectives as shown in Fig. 4b. It can be attributed to the fact that with the augment of decision variables, the situation of individuals not dominating each other gradually increases so that it is more necessary for such a method to accurately distinguish each individual. Furthermore, selecting the parent population according to the search bias can effectively decrease the loss of promising solutions and accelerate the convergence. Overall, the proposed strategy can markedly enhance the quality of the parent population for the expensive CMODOPs, especially when the problem size becomes larger.

Comparative experiment of model management strategy

To verify the effectiveness of the model management strategy used [54] on the proposed algorithm, two aspects of performance comparison are designed in this section. One is the RMSE between the predicted and the true values of an objective during optimization, as shown in Figs. 5 and 6, which is used to measure the evaluation accuracy of the RF model with or without the model management strategy, separately. The other is that the proposed algorithm is compared with RFMOISR without the model management strategy (RFMOISR-WMM) on the comprehensive performance including convergence and distribution, as shown in Table 2. Experiments were performed in MOKP instances with 30–50 decision variables and 2–3 objectives.

Fig. 5
figure 5

RMSE of the RF model in RFMOISR with or without model management strategy on the MOKP instances with 2 objectives

Full size image
Fig. 6
figure 6

RMSE of the RF model in RFMOISR with or without model management strategy on the MOKP instances with 3 objectives

Full size image
Table 2 Statistical results of the HV and IGD values on the two MOKP instances
Full size table

From Fig. 5, we can observe that, guided by the model management strategy, RMSE is gradually reduced to an acceptable error range during the iteration. However, the evaluation error of the RF model is always large without the use of model management, and even has an upward trend with the iteration. The above phenomenon clearly indicates that the model management strategy used has an excellent improvement effect on the prediction accuracy of the RF model in the proposed algorithm. Meanwhile, a similar conclusion can be drawn from the MOKP instance with 30–50 decision variables and 3 objectives in Fig. 6. This again illustrates that the model management strategy plays a crucial role in RFMOISR algorithm for solving CMODOPs.

In addition, the HV and IGD values of RFMOISR with or without the model management strategy in the four instances above are presented in Table 2. It can be seen from Table 2 that the proposed RFMOISR is significantly better than that without the model management strategy in terms of HV index and IGD index, which means that the model management strategy can bring outstanding convergence and diversity to the proposed algorithm. This can be attributed to the improvement of the evaluation accuracy of the surrogates by the model management, so that the search of the algorithm can be accurately guided. Moreover, coupled with the ISR strategy, the algorithm accelerates the convergence while maintaining the diversity of the population, resulting in a set of reliable solutions.

Comparative experiment of RFMOISR

Comparative analysis of IGD

The mean IGD values of the five different algorithms on small and large-scale MOKP instances are listed in Table 3. The best values among all the compared algorithms for each instance are marked in boldface. To judge whether there is a statistical significance between RFMOISR and other compared algorithms, the Wilcoxon rank sum test with significance level p = 0.05 is employed in Table 3, where symbols ‘+’, ‘−’, ‘=’ indicate that RFMOISR is significantly better than, significantly worse than, or comparable to the compared algorithms, respectively. Figures 7 and 8 draw the change graphs of IGD values obtained by the overall algorithms based on small and large MOKPs with 2 and 3 objectives, separately. According to the Wilcoxon rank sum test, the RFMOISR algorithm outperforms the others on most of the small-scale test problems, and achieves the optimal IGD values on all the large-scale MOKP instances. Moreover, we can also observe from Fig. 7 that RFCMOCO is the second best-performing algorithm among the five compared algorithms. It can be found that better performance on expensive MOKP instances is acquired from SAEAs, i.e., RFMOISR, RFCMOCO. Also, as shown in Fig. 8, despite the increasing scale of problems, these two SAEAs still maintain their performance advantages. The phenomenon can be well interpreted, which also shows the significance of SAEAs. Owing to costly and limited exact evaluations, the surrogate model with low cost and high precision undertakes the main evaluation task, which can save computing resources while obtaining higher quality solutions, especially on more complex large-scale problems.

Table 3 Statistical results of the IGD values on the small and large MOKP instances
Full size table
Fig. 7
figure 7

The IGD values for all the algorithms on the small-scale MOKPs with two or three objectives

Full size image
Fig. 8
figure 8

The IGD values for all the algorithms on the large-scale MOKPs with two or three objectives

Full size image

What is more, it is clear from Fig. 7 to 8 that SPEA2 is the best-performing algorithm among the three MOEAs in such instances. And NSGA-II performs very competitively to SPEA2. For MOEA/D, it is poor than that of its competitors, and it does not obtain an acceptable solution for any test problem. This may be since the PFs of MOKP test problems are irregular and discontinued, while a set of well-distributed weight vectors generated by MOEA/D cannot guarantee a good distribution of obtained solutions. By contrast, unlike MOEA/D where the search directions are fixed by weight vectors, SPEA2 can effectively avoid the search biased to a certain area in CMODOPs by introducing fitness assignment. It is also one of the main reasons why RFMOISR is the most effective algorithm in expensive MOKPs.

Comparative analysis of HV

For completeness, all algorithms have also been compared using the HV index. Table 3 presents the average HV results on such test problems, including the Wilcoxon rank sum test with significance level p = 0.05. Similar conclusions can be drawn from the results given in Table 4. These results in Table 4 indicate that RFMOISR can provide significantly better results than the other algorithms on an overwhelming majority of the small and large tested problems in light of the Wilcoxon’s test. Besides, to intuitively analyze the performance of the compared algorithms, the HV values for such algorithms on MOKPs with two and three objectives are severally plotted in Figs. 9 and 10. It is also observed from figures that RFMOISR achieves significantly better comprehensive performance on most of the test problems.

Table 4 Statistical results of the HV values on the small and large MOKP instances
Full size table
Fig. 9
figure 9

The HV values for all the algorithms on the small-scale MOKPs with two or three objectives

Full size image
Fig. 10
figure 10

The HV values for all the algorithms on the large-scale MOKPs with two or three objectives

Full size image

There are two main reasons why the proposed algorithm shows comparably well on the HV index. One is the usage of RF models, and the other is the improved stochastic ranking approach. Random forest is adopted as the surrogate model according to its high approximation accuracy for discrete problems, which can reduce the evaluation error of individuals and correctly grasp the search direction. Moreover, the fitness mechanism takes into account both the dominant relationship and the aggregate distribution of population to accurately identifying the identify the potential of each individual, thereby improving the quality of parent individuals. And adaptive probability operator is developed in ISR to be consistent with the search principles of evolutionary algorithms to speed up the search process and converge quickly to an optimum. Hence, thanks to the novel mechanisms proposed, the RFMOISR algorithm still has significant potential for the large-scale MOKP instances in Fig. 10. Overall, the statistical analyses of the above two indexes have evidenced that RFMOISR can better balance convergence and diversity in expensive CMODOPs.

Comparative analysis of ME

In this subsection, the convergence of the compared algorithms can be observed from Table 5 which shows the mean and variance of the convergence metric ME obtained using five algorithms RFMOISR, RFCMOCO, SPEA2, NSGA-II, and MOEA/D. The Wilcoxon test is also adopted to judge the significance of the results compared to RFMOISR algorithm. The results of Table 5 show that RFMOISR is able to converge better on such MOKP test problems, which indicate SAEAs with the help of random forest models can refrain from a misleading search on expensive CMODOPs. Also, the collaboration of fitness assignment and adaptive selection can further contribute to improving the quality of solution and convergence speed. Furthermore, with the increase of problem scale and solving difficulty, the superiority of the proposed strategies under limited real evaluations are ulteriorly revealed, so that the algorithm can converge to a satisfactory level. Additionally, in all cases with RFMOISR, the variance in 30 runs is also small, which manifests the stability of the proposed algorithm in optimizing expensive constrained multi-objective discrete optimization problems.

Table 5 Statistical results of the ME values on the small and large MOKP instances
Full size table

Comparative analysis of GD

Since the algorithms failed to converge to the Pareto fronts on some of the benchmarks due to the limited evaluation budget, GD metric is adopted to investigate the convergence of the five algorithms. The statistical results in terms of GD values obtained by the five algorithms on small and large-scale MOKPs are summarized in Table 6. The best result on each instance is highlighted in bold. Also, the Wilcoxon rank sum test is also employed to compare the results at a significance level of 0.05. Symbol +, −, and ≈ means the result of RFMOISR is significantly better than, worse than or similar to that of the compared algorithm as well.

Table 6 Statistical results of the GD values on the small and large-scale MOKP instances
Full size table

As listed in Table 6, the proposed RFMOISR obtains the best GD values compared to other four algorithms on a vast majority of the small-scale test problems. And it achieves significant advantages with respect to the convergence in large-scale MOKP problems, confirming the better performance indicated by ME values. Overall, it can be concluded that the proposed algorithm has the preeminent capability in approximating Pareto fronts of such expensive CMODOPs.

Comparative analysis of computational cost

In this section, we perform a comparative analysis of the computational cost of RFMOISR and RFCMOCO to gain some insight into the computation efficiency of the proposed method. Table 7 presents the average computation time consumed by the competitive algorithms over 30 independent runs on the ten small-scale and six large-scale benchmark problems with 10–750 decision variables under the limited computational budget. From Table 7 together with Table 3, we can find that the proposed RFMOISR requires less computation time whilst obtaining better results than the compared algorithm RFCMOCO in 6 out of 10 small-scale cases. Moreover, the RFMOISR outperforms the RFCMOCO on the large-scale test problems in terms of solution performance and running time. This means that the ISR strategy, while improving the quality of the solutions, can also notably save the computational overhead of RFMOISR and accelerate the search efficiency of the algorithm.

Table 7 Statistical results of the runtime(s) values on the MOKP instances
Full size table

For visual observations, Figs. 11 and 12 draw the runtime of the two SAEAs for the small- and large-scale MOKPs with two or three objectives, respectively. The graphs show that RFMOISR is significantly more efficient than RFCMOCO in such sixteen benchmark problems. Specifically, it can be seen from Fig. 11a that compared with RFCMOCO, the proposed algorithm reduces the calculation time by at least 32% on the test functions with two objectives. What is more, Fig. 11b manifests that the proposed algorithm can be saved at least 50% of the calculation time of RFCMOCO in such MOKP instances with three objectives. This again reinforces the advantages of the improved stochastic ranking, which has the nature of computation simply and provides a good parent population for the next generation, thereby further saving computing resources. As it saves some evaluation costs, it gets more opportunity to explore the search space. Although the size of the problem becomes large in Fig. 12, the proposed algorithm still performs prominently in convergence speed due to the above strategies. To sum up, it is potential for RFMOISR to obtain better solutions in a relatively short time on expensive constrained multi-objective discrete optimization problems.

Fig. 11
figure 11

The runtime of the two SAEAs on the small MOKPs with 2 or 3 objectives

Full size image
Fig. 12
figure 12

The runtime of the two SAEAs on the large MOKPs with 2 or 3 objectives

Full size image

Besides, we find from Fig. 13 that as the number of objectives and decision variables increases, the complexity of the problem also rises, so the time taken by the two SAEAs augments apparently. However, the time taken by RFMOISR to optimize the instances with three objectives is close to or even less than the time taken on the two objectives when the problem scale is increased to more than 30. This may illustrate that the proposed algorithm is more efficient in solving the expensive high-dimensional CMODOPs.

Fig. 13
figure 13

The runtime of the two SAEAs on all the MOKPs with 2 and 3 objectives

Full size image

Conclusion and future work

In this paper, we propose a surrogate assisted multi-objective optimization with improved stochastic ranking, RFMOISR, for computationally expensive constrained multi-objective discrete optimization problems, in which the random forest surrogate models are coordinated with an individual-based model management strategy to enhance approximation accuracy. Moreover, an improved stochastic ranking method is proposed by introducing the fitness mechanism and adaptive probability operator to balance better objectives and feasible solutions and improve the quality of the parent population for the next generation, which is a simple and efficient approach to compensate for the time-consuming nature of the RF model. Experimental results on the ten small-scale and six large-scale typical benchmark problems demonstrate the performance of the proposed RFMOISR algorithm for expensive constrained multi-objective discrete optimization problems outperforms the existing schemes in terms of convergence, diversity and computational cost.

In the experimental analysis, the improved stochastic ranking strategy is compared with original stochastic ranking, stochastic ranking with adaptive selection, and stochastic ranking with the fitness mechanism to verify the superiority of its selection accuracy based on the RF model, especially on expensive constrained multi-objective discrete optimization problems with larger scales. Meanwhile, the significance of individual-based model management to the proposed algorithm is verified from two aspects of approximate error and algorithm performance. Furthermore, different kinds of performance metrics were employed to check the effectiveness and quality of the proposed algorithm. The statistical results obtained from RFMOISR when compared to other popular algorithms including RFCMOCO [38], NSGA-II [55], SPEA2 [53], and MOEA/D [56] indicated that solutions of RFMOISR are with the highest quality in terms of both the convergence and distribution under the limited computing resources. Besides, compared with the RFCMOCO in running time index, we also find that the proposed algorithm can save more computational cost. To sum up, the proposed algorithm is promising to address the expensive constrained multi-objective discrete optimization problems.

In the future, we are planning to investigate more efficient approaches for solving expensive constrained multi-objective combinatorial optimization problems. Also, RFMOISR is expected to be applied to more complicated engineering problems to further prove its potential for real-world applications.