A comparison among optimization software to solve bi-objective sectorization problem

In this study, we compare the performance of optimization software to solve the bi-objective sectorization problem. The used solution method is based on an approach that has not been used before in the literature on sectorization, in which, the bi-objective model is transformed into single-objective ones, whose results are regarded as ideal points for the objective functions in the bi-objective model. Anti-ideal points are also searched similarly. Then, using the ideal and anti-ideal points, the bi-objective model is redefined as a single-objective one and solved. The difficulties of solving the models, which are basically non-linear, are discussed. Furthermore, the models are linearized, in which case how the number of variables and constraints changes is discussed. Mathematical models are implemented in Python's Pulp library, Lingo, IBM ILOG CPLEX Optimization Studio, and GAMS software, and the obtained results are presented. Furthermore, metaheuristics available in Python's Pymoo library are utilized to solve the models' single- and bi-objective versions. In the experimental results section, benchmarks of different sizes are derived for the problem, and the results are presented. It is observed that the solvers do not perform satisfactorily in solving models; of all of them, GAMS achieves the best results. The utilized metaheuristics from the Pymoo library gain feasible results in reasonable times. In the conclusion section, suggestions are given for solving similar problems. Furthermore, this article summarizes the managerial applications of the sectorization problems.


Introduction
Almost all optimization software and solvers have evolved considerably over time, but they still have significant differences [1]. Therefore, in the literature, there are studies comparing solvers [2][3][4][5][6][7]. Likewise, some studies compare the evolutionary algorithm (EA) with solvers and discuss the advantages and disadvantages [8][9][10]. In this study, a comparison is made between optimization software to solve a sectorization problem (SP). The goal of SPs is to split a large region into smaller ones for administrative purposes [11][12][13][14]. However the main contribution of the study to the literature is the comparison among optimization software, it has other novelties, which can be summarized as follows: • A new bi-objective (BO) model is defined for SP, which is nonlinear (NL). • A new approach is applied to solve the model, in which, at first, the BO model is divided into two single-objective (SO) ones, which are also NL. • The goal of the division is to gain ideal points (IDPs) and anti-ideal points (ADPs) for the objective functions in the BO model.  • The results of the linear (L) and NL models are supposed to be the same if they are optimal. However, this is also a way to evaluate the computational performance of the solvers. • The BO model is re-defined based on the LP-metric method, using IDPs and ADPs for the objective functions.
Though the employed solution approach is well-known in optimization topics, it is employed for the first time in the sectorization literature.
As a contribution to the literature, a comparison is made between different software for solving MO SP. For mathematical models, different solvers like Python's Pulp library, Lingo, IBM ILOG CPLEX Optimization Studio, and GAMS are used. There is a similar study in the literature comparing optimization tools for solving SPs [15]. But it doesn't report the results of GAMS, IBM ILOG CPLEX Optimization Studio, and Lingo. Furthermore, the solution approach in this study is different. From SO metaheuristics available in the Pymoo library of Python, a Genetic Algorithm (GA), and from MO ones the Non-dominated Sorting Algorithm (NSGA-II) are employed [16].
A summary of the methods that are applied in this paper is presented in Fig. 1, whose details are described in Section 4. The used abbreviations can be found in the appendix.
The other parts of the paper are organized as follows: Section 2 summarizes the related literature. The problem definition and solution method are respectively described in Sections 3 and 4. In the experimental results section, based on generated benchmarks, the operated software and algorithms are compared. Conclusion and future works form the last part of the article, where comments and suggestions for solving similar models are offered as well as several research directions.

Literature review
Koch et al. [1] report that from 2001 to 2020 linear programming (LP) and mixed integer linear programming (MILP) solvers became about 180 and 1000 times faster, respectively. Furthermore, benchmarks that couldn't be solved years ago are now solved in seconds [1]. Despite all these advancements, occasionally unexpected variability affects the performance of mixed integer programming (MIP) solvers. This issue has been observed for decades. Lodi and Tramontani [2] discuss the causes of this phenomenon. Therefore, the performance of solvers for various problems is of interest to researchers. Neumaier et al. [3] compare the performance of constrained optimization solvers over 1000 benchmarks from literature, which contain 1000 variables. Kronqvist et al. [4] make a review and comparison between mixed integer non-linear programming (MINLP) solvers based on a set of 335 benchmarks, which are divided into subsets based on the continuous relaxation gap, the degree of nonlinearity, and also the relative number of discrete. Table 1 Some of the problems modeled based on sectorization and employed approaches to solve them.

Problem
Solution Approach MO model Political districting [11] I L P Facility location [12] G P ✓ Location allocation [13] Metaheuristics Location routing [17] Heuristic and metaheuristic Location routing [18] H e u r i s t i c Supply chain network [19] Metaheuristic Elderly care service districting [20] M I N L P ✓ Emergency medical service districting [21] Metaheuristic Power distribution networks [22] Metaheuristic ✓ Police districting [23] Literature review Commercial districting [24] Metaheuristic School districting [25] H e u r i s t i c Air traffic control [26] Heuristic and metaheuristic Air traffic control [27] Human factors issues Air traffic control [28] Constraint programming Air traffic control [29] Metaheuristic ✓ Air traffic control [30] Simulation Air traffic control [31] Metaheuristic ✓ Telecommunications [32] Simulation Telecommunications [33] I Q P , M I P Water distribution networks [34] MINLP and metaheuristic ✓ Table 2 Managerial applications of SPs.

Application Filed
Decreasing wastes [35] Water management Managing capacity [36] Air traffic management Improving utility [37] Healthcare management The acquired results provide information about which solver or method is suitable for which type of benchmarks. Meindl and Templ [5] compare various solvers on LP benchmarks. Luppold et al. [6] investigate the effects of different types of variables, constraints, and logical expressions on the performance of solvers. Describing the complexity of hard real-time systems with large numbers of constraints, Guasque and Balbastre [7] emphasize the importance of heuristic methods to obtain feasible solutions for them. They also make a comparison between CPLEX and Gurobi. In this study, the comparison between optimization software is done over SPs. As stated in the previous paragraph, the type of problem and the solution method affect the performance of employed solvers. Therefore, this section summarizes the literature on different types of SPs and the utilized methods for their solution. The main objective of SPs is to divide a vast district into smaller and balanced areas. SPs are mostly formulated as multi-objective (MO) models and they have applications in numerous fields. Some of the problems modelled based on sectorization, the used solution approaches for them, and whether the model is MO or not appear in Table 1. The used abbreviations can be found in the appendix.
SPs have managerial implications, some of which are shown in Table 2. MO models fall into multi-criteria decision-making (MCDM) techniques. As seen in Fig. 2, MCDM techniques can be categorized as multi-attribute decision-making (MADM) and multi-objective decision-making (MODM) subcategories that respectively have discrete and continuous decision spaces. MODMs are divided into four categories based on the role of the decision maker (DM). In the first class, which is a neutral approach, there is no preference information from the decision-maker. The LP-metric method is in this category, which trade-offs between desired and undesired solutions and in this way converts an MO model to an SO one. In other classes, preference information is received from the DM before, after, or interactively during the resolution process. MADM techniques can be divided into three categories: compensatory, non-compensatory, and partial compensatory. In compensatory approaches, some negative attributes can be compensated by positive ones [38][39][40][41].
In the literature, the LP-metric or the methods designed based on IDPs and ADPs have not been used to solve SPs, however, they have been applied for the solution of several problems, some of which are summarized in Table 3 [42].
Some of the problems in the literature solved based on other techniques of MCDM are summarized in Table 4. There are also other MCDM techniques in the literature that have not been used for SPs. For example, from MODM, a priori methods: lexicographic max-min [53] -constraint, goal programming [54], surrogate worth trade-off method (SWTM) [55], from interactive methods [56], step method (STEM) [57], from MADM, compensatory methods: the technique for order preference by similarity to ideal solution (TOPSIS), analytical network process (ANP) [58][59][60][61] and from MADM, non-compensatory methods: Max-Max, Min-Max [62] have not been employed for SPs solution.

Table 3
Problems solved based on IDPs and ADPs in the literature.

Problem
Combined Method with IDPs and ADPs Location selection [43] Fuzzy decision making Efficiency evaluation [44] Data envelopment analysis Supplier selection [45,46] M C D M Wind site selection [47] Fuzzy MCDM Portfolio asset management [48] MODM Table 4 Some of the problems solved based on MCDM techniques.

Problem Method
Category of the Method Air traffic control [49] Weighted sum MODM, a priori methods Water distribution network [50] Lexicographic MODM, a priori methods Air traffic control [51] Scalarization MODM, a priori methods Water distribution network [52] AHP MADM, compensatory

Problem definition
In the problem, it is assumed that there are clients, to be assigned to some service centers, whose numbers and coordinates are predetermined. Each service center and its assigned clients form a sector. Each client is assigned to only one sector. Each of the formed sectors has at least one client. The number, coordinates, and service centers are known beforehand. Each client has a predefined demand.
Some of the used notations are summarized in Table 5.
The binary decision variable of the model is as in Eq. (1).

BO model
To be minimized, the objective function related to the equilibrium of distances in sectors is defined as in Eq. (2).
. To be minimized, the objective function related to equilibrium of demands in sectors is defined as in Eq. (3).
. The purpose of 1 and 2 is to create more balanced sectors in terms of distance and demand. In this way, to be minimized, the BO function of the model is defined as in Eq. (4).
Each client is assigned to just one sector, which is satisfied by Constraint (5).
As specified in Constraint (6), at least one client is assigned per sector.
Furthermore, as defined in Eq. (1), is a binary variable, which is considered as a constraint in the model. It is expected that the equilibrium is yielded in terms of distances and demands between the formed sectors, which is provided by Constraints (7)- (8).
, and are parameters for the equilibrium of distances and demand in sectors. The reason for using Constraints (7)-(8) can be explained as follows: The total deviations of and , = 1, ..., , from their averages are minimized by applying the objective functions 1 , and 2 . Whereas using Constraints (7)-(8), more desired solutions are provided by defining upper limits to the deviation value for each sector. If the objective functions were defined based on variance or standard deviation instead of absolute value, they would have the same influence, but upper limits of deviations could not be guaranteed.
The complete form of the basic NL model can be expressed with Statements (9)- (15).

Solution method
Before describing the solution method, some concepts are explained. Vectors of nadir [63] and ideal objectives are used in the literature on MO optimization [64], which are not the same as IDPs and ADPs in this work. The next definitions need to be made to clarify the differences: In a minimization problem the following conditions must be valid for the solution 1 to dominate 2 :

Definition 4.2.
is a Pareto optimal or an efficient solution if it is not dominated by any other solution.
In the field of MO optimization, two distinct terms, "weak efficiency" and "strict efficiency," are distinguished. A solution is considered weakly efficient if it is not dominated by another one in the solution space. This means that, given a certain objective, a weakly efficient solution may be dominated by another one. Yet, a solution is considered strictly efficient if it is not dominated by any other solution according to any objective. In this research, weakly efficient solutions are referred to as efficient solutions. The IDPs and ADPs in this study are the outputs obtained from SO models, which are operated in the LP-metric method to transform the MO model into an SO one.
It can be stated that the LP-metric method sets weights for objectives and hence is similar to the weighted sum method. The weighted sum approach is in the category of a priori methods in Fig. 2, because weights are usually defined by the DM before the solution process. Whilst in the LP-metric method, the weights are determined during the solution procedure.
According to Ref. [38], optimal solutions of the weighted sum problem with positive (non-negative) weights are always (weakly) efficient. Furthermore, all (weakly) efficient solutions are optimal solutions of scalarized problems with positive (non-negative) weights under convexity assumptions.
It should be noted that integer programming models are considered non-convex, thus convexity assumptions are not valid for them.
In the next subsection, the solution method is described, which has not been used in previous studies about SP. In this method, at first, the BO model is converted into two SO models, from which IDPs and ADPs are found for the objective functions. Then the BO model is redefined and transformed to an SO model by using IDPs and ADPs. L versions of the models are also utilized.

Dividing the BO model into two NL SO models
The goal of the division is to determine IDPs and ADPs for objective functions.

4.1.1.
1 , SO model to find an IDP for 1 in the BO model To be minimized, the objective function of this model is 1 , as defined in Eq. (2). The optimal solution of this model, * 1 , is an IDP for 1 in the BO model.

4.1.2.
2 , SO model to find an IDP for 2 in the BO model To be minimized, the objective function of this model is 2 , as defined in Eq. (3). The optimal solution of this model, * 2 , is an IDP for 2 in the BO model.

Linearizing SO models
In this subsection, the SO models are linearized according to two different approaches. The number of variables in these two approaches is different. They are used to demonstrate the influence of the number of variables.

4.2.1.
1 , linearized SO model based on the first approach to find an IDP for 1 in the BO model Using non-negative variable 1 objective function 1 is linearized.

4.2.2.
1 , linearized SO model with the second approach to find an IDP for 1 in the BO model Using non-negative variables 1 and 2 , = 1, ..., , objective function 1 is linearized.
To be minimized, the objective function of this model is determined as in Eq. (24). The optimal solution of this model, denoted as * 1 , is an IDP for 1 in the BO model. Constraints (5)- (6) and (20)- (23) are also contained in this model, as well as Constraints (25)- (26).
In the 1 model, there are additionally 2 variables due to 1 and 2 , while there are further variables in the 1 model because of 1 or 2 . Moreover, in both models, 7 constraints are added.

4.2.3.
2 , linearized SO model based on the first approach to find an IDP for 2 in the BO model Using non-negative variable 2 objective function 2 is linearized.
Constraints (5)-(6), and (20)- (23) are also included in this model. The optimal solution of this model indicated as * 2 , is an IDP for 2 in the BO model.

4.2.4.
2 , linearized SO model based on the second approach to find an IDP for 2 in the BO model Using non-negative variables 3 and 4 , = 1, ..., , objective function 2 is linearized. To be minimized, the objective function of this model is determined as in Eq. (31), taking Constraints (32)-(33) into consideration, as well as Constraints (5)-(6), (20)-(23). The optimal solution of this model denoted as * 2 , is an IDP for 2 in the BO model.
The number of added variables and constraints in the 2 and 2 models is the same as described at the end of Subsection 4.2.2.

4.2.5.
1 , SO model with linearized objective function based on the first approach to find an ADP for 1 in the BO model The objective function of this model is defined as in Eq. (16), to be maximized, considering Constraints (17)- (19), in addition to Constraints (5)- (6). Specified as * * 1 , the optimal solution of this model is an ADP for 1 in the BO model.  (30). Shown as * * 2 , the optimal solution of this model is an ADP for 2 in the BO model.

4.2.8.
2 , SO model with linearized objective function based on the second approach to find an ADP for 2 in the BO model To be maximized, the objective function of this model is specified as in Eq. (31), regarding Constraints (32)- (33), in addition to Constraints (5)- (6). Indicated as * * 2 , the optimal solution of this model is an ADP for 2 in the BO model. As mentioned before, since Constraints (7)-(8) identify upper limits, are not applicable for maximization. Therefore, the and models don't include them.

Redefining the BO model
Since it is not easy to deal with multiple objective functions simultaneously, in this section the BO model is redefined as an SO model using IDPs and ADPs.
The LP-metric method is used to convert an MO problem to an SO one, in which the total distance to the IDPs is minimized, along with maximizing the distance from the ADPs.
In the next subsection, three models are defined based on the LP-metric method.

LP-metric-NL, an NL model using the outputs of the and models and LP-metric method to find efficient solutions for the BO model
Using the results obtained from the NL models, to be minimized, the objective function can be defined as in Eq. (34).
where 1 and 2 are as in Eqs. (2) and (3) and * , * * are the outputs of the models and , = 1, 2, respectively. In this model Constraints (5)-(8) are regarded. Applying the outcomes achieved from the L models based on the first approach of linearization, to be minimized, the objective function can be determined as in Eq. (35).

LP-metric-, a L model using the outputs of the and models and LP-metric method to find efficient solutions for the BO model
Using the results acquired from the L models based on the second approach of linearization, to be minimized, the objective function can be determined as in Eq. (35), in which 1 and 2 are as in Eqs. (24) and (31), and * , * * are the outputs of the models and , = 1, 2, respectively. Constraints (5)-(6), (20)-(23), (25)-(26), (32)- (33) are considered in this model. For = 1, = 2, and = ∞, the distances are the street-block (or Manhattan), Euclidean, and Chebyshev. A Pareto optimal solution is obtained for each . But the higher , the higher the complexity. In this study, only = 1 is used to simplify.

Experimental results
We generate four benchmarks as 200×10, 500×31, 1000×76 and 2000×200, which are pointed out as the × . = ⌊ ⌋ is the target number of points per sector. In order to have diversity in the benchmarks, Q is equal to 20, 16, 13, and 10, respectively. Parameters ranges inside the interval [0, 1]. Two-dimensional coordinates of clients and center of sectors and likewise clients' demands, are created according to (50;10) and (10;100), which are normal and discrete uniform distributions, respectively [13,14,17,18].
The proposed models are not very complex, but the solution becomes difficult because the models are basically NL and also the number of variables and constraints in the benchmarks can be high, especially in the L models. For example, in the NL models, the number of the binary variables is × . Because of the Constraints (5), (6), (7), (8) respectively, linear, linear, non-linear and non-linear constraints are included in the model. As an example, the total number of variables and constraints for the model 1 are given in Table 6, where and stand for the total number of variables and constraints, respectively. Subsections 4.2.2 and 4.2.4 specify how the number of variables and constraints varies in the linear models. A system with an Intel Core i5 processor, 2.4 GHz with 12 GB of RAM is utilized. IBM ILOG CPLEX Optimization Studio 12.9, GAMS 25.1.2, Lingo 19.0, Pymoo 4.2, and Pulp 2.4 are operated. From GAMS solvers BARON is used for the NL models, while CPLEX is used for L ones. Even if BARON is used for L models it is shifted to CPLEX automatically. In the tables of results, CPLEX refers to IBM ILOG CPLEX Optimization Studio, not the solver in GAMS. For the GA and NSGA-II utilized from the Pymoo library, we set the size of the population and offspring equal to 1000. We use a random sampling method and uniform crossover and mutation [16].
We suppose 0 ≤ ≤ 1, and 0 ≤ ≤ 1. Using Constraints (7)- (8), the values of and can be calculated as in Statements (36)-(37): We employ GA to define the values of these parameters. Such that, running the GA for the model SI, starting from 1 and decreasing by 0.05 each time, the highest value that gives a feasible result is selected. We assign the values of these two parameters equally to manage them easily. Thus, the acquired values are as in Table 7. The reason for using GA for this purpose is that it provides results faster.
The results are presented in Tables 8-14. Using the same distributions we generate alternative instances to verify whether the acquired results depend on the employed benchmarks. Obtained outputs for the alternative benchmarks are given in Table 10. But the results in Tables 8-9 and 12-14 are for the main benchmarks.  Although the models are not complex, it is not easy to find their optimal solutions. This is because they are basically NL and when they are linearized the number of variables and constraints grows. Therefore, the tables of results mostly report a feasible solution.

A. Teymourifar
With GA and NSGA-II metaheuristics, only NL models are solved, because the linearization is done to reduce the complexity of the problem, which is not a handicap for these algorithms and they can easily solve the NL models. In the tables, stands for solution time in a format as minutes:seconds (mm:ss). Optimal results are displayed as . The maximum working time of the solvers is defined as 10 minutes. If no solution is obtained within this time, it is shown in the tables as NaN. If only a feasible solution is found in 10 minutes, this is signed by .
As seen in Table 8, only GA is able to find solutions for all benchmarks of the model SI within the defined time, though their quality is not satisfactory. Although Pulp and Lingo are capable of solving L models, they only find results for the benchmark with the smallest size. Hence, their outcomes are not presented in Table 9. In the L models, although GAMS uses the CPLEX solver, as seen in Tables 9-10, there is a difference between its results and the IBM ILOG CPLEX Optimization Studio. It is observed that GAMS achieves better results for more benchmarks. Moreover, according to the results, linearization based on the first and second approaches does not seem to predominate each other. In Tables 8-10 the best achieved solution for each model is bold.
To understand if there are significant differences between the results of the solvers, we perform an analysis of variance (ANOVA) with a confidence level of 95%. We use the outputs of Table 9. We define factors as the type of solver, size of the benchmark, and model. So, the analysis is a three-way ANOVA. The null hypothesis is that the results of groups of factors are equal, against the alternative hypothesis that at least one of them is different.
As seen in Table 11, since the p-values for the factors of solver and size are less than 0.05, they are significant. However, the p-value for the factor of the model is greater than 0.05, and consequently, it is not significant. This is also valid for the interactions of the factor of the model with other ones, which are shown as Solver:Model and Size:Model.
While the ANOVA indicates a significant overall difference, it does not determine which specific groups differ from each other. To compare groups, we perform a multiple comparison test, which is a post hoc analysis. The p-value for the test is equal to zero. Two groups are formed based on the results of GAMS and IBM ILOG CPLEX Optimization Studio solvers, whose marginal means of groups are significantly different. Overall, based on the results of the ANOVA and multiple comparison tests, we can conclude that GAMS and IBM ILOG CPLEX Optimization Studio perform differently.
As seen in Fig. 3a-d, except for Model 1 of benchmark 500 × 31, IBM ILOG CPLEX has better performance than GAMS. But as it appears in Table 10, GAMS outperforms IBM ILOG CPLEX for other benchmarks.  Table 9 Results obtained for the model LSI.  A. Teymourifar Table 11 The results of the three-way ANOVA for the outputs of Table 9.

Table 14
Obtained solutions for the benchmarks of the MO model using LP-metric-NL, LP-metric-and LP-metric-methods.
Benchmark LP-metric-NL LP-metric-LP-metric- It should be noted that for the LP relaxation of the LSI model, in all benchmarks, both GAMS and IBM ILOG CPLEX Optimization Studio find the optimal result in less than one minute. What is meant by LP relaxation is just that the binary decision variable in Eq. (1) is defined continuously in the interval [0, 1].
GAMS, IBM ILOG CPLEX Optimization Studio, and Lingo could not provide solutions for the SA and LSA models, reporting that the models are unbound. It may be useful to apply additional constraints to avoid this, but this is not the subject of this study. Because the aim of this work is to solve the presented model with existing constraints. Nevertheless, the employed GA gives results for SA as in Table 12, which are supposed as ADPs.
As mentioned before, the NL types of the BO model are also solved with NSGA-II available in the Pymoo library of Python. The results are given in Table 13, where shows the number of obtained non-dominated solutions. It has been mentioned earlier that the optimal solutions of the SO models will be utilized as IDPs and ADPs for the BO model, while according to the acquired results, most of them are from the feasible ones. This can be regarded as a limit for the study. Using the IDPs ve ADPs in the LP-metric method, defining = 1, outputs are gained as in Table 14.
In addition to 1 and 2 , the values of are also given, which is defined as in Eqs. (34) and (35). As it has been mentioned earlier, the results in Table 12 are regarded as ADPs in all models and * * = * * , = 1, 2. The LP-metric-NL is solved with GA while the LP-metric-and LP-metric-are handled with GAMS. Related outputs are presented in Table 14 and Fig. 4, where non-dominated solutions are bold. As seen, some of the non-dominated solutions achieved by the LP-metric method for benchmarks 200×10, 500×31, 1000×76 are the same as those in Table 12. Even for benchmarks 2000×200, the solution provided by the LPmetric-NL method dominates the ones in Table 12. Though, the results of the LP-metric-and LP-metric-are significantly better than the LP-metric-NL, which is also due to the performance of the CPLEX solver in GAMS.
A. Teymourifar   Fig. 4. Visualization of the results presented in Table 14. As mentioned earlier, the reason for the definition as = 1 is to don't ascend the complexity of the problem. But consequently, some of the solutions on the Pareto front are overlooked. It is possible to find other non-dominated solutions with different values of [65]. The models designed based on the LP-metric method can only achieve a few of the supported non-dominant solutions. Such solutions are located on the smallest convex polyhedron that includes the feasible region, which is called the convex hull. Whereas the unsupported ones are inside the convex hull. Examples of these solutions are depicted in Fig. 5. Unlike the unsupported nondominant solutions, the supported one can be uncovered by applying weighted sum or similar methods like LP-metric. Discovering unsupported non-dominant solutions is a challenge for mathematical programming techniques [66][67][68]. Based on the similarity between the LP-metric approach and the weighted sum method, it can be emphasized that it is not always possible to identify all non-dominated solutions using this approach.

Conclusion and future research directions
In this study, for the first time in the literature, a comparison is made between various solvers like Lingo, IBM ILOG CPLEX Optimization Studio, GAMS, and also Pulp, and Pymoo libraries from Python for solving SP. The comparison is made based on a new BO model for an SP. The difficulties of solving the model, which is basically NL, are discussed. The effects of objective functions and the non-linearity of some constraints on this difficulty are mentioned. An approach that had not been used in the literature before is used to solve the model, in which, at first, the BO model is divided into two NL SO models. The results of the SO models, when the objective function is to be minimized and maximized, are assumed to form IDPs and ADPs for the respective objective function in the BO model. Also, the models, which are basically NL, are linearized to solve the BO model easier. Both mathematical programming and metaheuristics are used for the solution. Although the models in the study are not complex, solutions are difficult due to nonlinearity and also because of the large number of variables in some benchmarks when they are linearized. For these reasons, solvers are not very efficient in solving models, and in general, during the defined run time, they could not find optimal solutions. Solvers performed better on the L models than the NL ones. Usually, GAMS is known to be one of the best software for operational research. In this study, GAMS provides better results than Lingo and Pulp. Even though GAMS operates the same solver as IBM ILOG CPLEX Optimization Studio for L models, it is slightly better. Among the solvers, the weakest performance is for Pulp.
The applied SO and MO metaheuristics are able to yield feasible solutions in reasonable times though they couldn't find optimal results. Their run times are much less than the other solvers but the quality of the results is not satisfactory. One of the reasons for this matter can be that the used metaheuristics are from a generic open-source library and they are not designed problem-oriented.  Nevertheless, it may be useful to use the algorithms of such libraries for solving large-scale problems, where it is not easy to find optimal solutions.
The models defined based on the LP-metric method achieve good results in terms of providing a non-dominated solution for each benchmark. Therefore, it would be said that the approach used for the first time in the literature of MO SP can be useful in finding non-dominant solutions. If more types of distance are used in the formulation, it could be possible to acquire more non-dominated solutions, but in this case, the complexity of the models would increase. This can be seen as a limit for this method. Notwithstanding, it can be concluded that the approach suggested in this study for solving MO SP is practical. The approach can be summarized as follows: dividing the BO model into some SOs ones, finding IDP and ADP for each objective function, reformulating and solving the BO model, and using the IDPs and ADPs in the LP-metric approach.
A limitation of this study is that IDPs and ADPs are got from different approaches. It is reported by the solvers that the maximization of the models is unbound. Although this matter is not an impediment to the completion of this work, it showed that it is better to reformulate the models, in future studies.
Figuring out the deficiencies of the metaheuristics employed in this study, in future studies, it is planned to propose a metaheuristic specific to solve MO SPs, in order to cope with the difficulties arising from a large number of variables and non-linearity in the models.
Performance metrics are used for quantifying performance and comparison of MO methods. It is crucial to use metrics to do comparisons between MO EAs designed based on the concept of Pareto optimality. In most MO problems, the Pareto optimal set is unattainable, correspondingly the measure of performance is challenging. Pareto front generation in MO problems needs high computational effort, especially in large-size problems. So its formation in a short time and with a small memory represents the efficiency of the employed algorithm. Convergence, coverage, and success rate measures are used to compare MO algorithms designed based on the Pareto approach. There are also other measures based on accuracy and variance. Accuracy shows generated non-dominated and the best-known solutions closeness. Coverage and variance have alike definitions. Coverage indicates solutions diversity and distribution. Variance represents the non-dominated front maximum range, which is covered by the generated solutions, for each objective. Diversity-based criteria measure the distances of the solutions in the first front of the final population. This front solution set is compared with a uniform distribution and its deviation is calculated. If an algorithm achieves more uniformly distributed solutions with less average deviation from the first front, it has good performance. Furthermore, criteria are categorized as quality, spacing, and diversification metrics. The obtained non-dominated solutions of algorithms are put together and their ratio is calculated as the quality metric. Spacing metric measures the uniformity of the solutions, while the diversification metric measures the spread of the solution set. There are also metrics such as net front contribution ratio and mean relative percentage increase to measure non-dominated solution set quality [69][70][71][72][73][74][75][76].
A significant limitation of the study is that the comparisons are over limited benchmarks and without a comprehensive statistical design of experiments [13]. In particular, this matter gains importance as the performance variability of the solvers is evident in the results. Likewise, in addition to SPs, benchmarks of similar problems such as p-median [77] should also be contained in the comparison. Statistical analyzes should be performed within a certain confidence interval, taking into account the different levels of various factors which can be problem type, benchmark size, solution, and linearization method [78]. The impacts of factors can be evaluated by factor analysis. The discrepancy between solvers' performance should be explored by designing hypotheses, performing ANOVA, and also post hoc methods such as Tukey's test [79,80]. A possible framework is summarized in Fig. 6.
Currently, there are not enough studies about the performance of solvers in the literature. In case of existing enough studies on this subject, a meta-analysis can be done, which is a research approach that can explore publications' bias and can draw a broader conclusion. It employs statistical methods to quantitatively compare and also combines the results of multiple studies [81,82].
In the literature, not enough attention has been given to the MADM techniques outlined in Fig. 2 to solve SPs. In real-life problems, DMs usually assign dissimilar priorities to criteria, in which case both compensatory and non-compensatory MADM can be beneficial.
Integrating MIP techniques with neural networks can provide an intelligence search on the feasible space [83,84], which will be considered in future works.

CRediT authorship contribution statement
Aydin Teymourifar: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement
Data will be made available on request.

Acknowledgements
This work was supported by The Portuguese Foundation for Science and Technology (FCT) in the framework of the project CEECINST/00137/2018.
The financial support from FCT (through the project UIDB/00731/2020) is also gratefully acknowledged.