Research Article  Open Access
Laurent Lemarchand, Damien Massé, Pascal Rebreyend, Johan Håkansson, "Multiobjective Optimization for Multimode Transportation Problems", Advances in Operations Research, vol. 2018, Article ID 8720643, 13 pages, 2018. https://doi.org/10.1155/2018/8720643
Multiobjective Optimization for Multimode Transportation Problems
Abstract
We propose modelling for a facilities localization problem in the context of multimode transportation. The applicative goal is to locate service facilities such as schools or hospitals while optimizing the different transportation modes to these facilities. We formalize the School Problem and solve it first exactly using an adapted constraint multiobjective method. Because of the size of the instances considered, we have also explored the use of heuristic methods based on evolutionary multiobjective frameworks, namely, NSGA2 and a modified version of PAES. Those methods are mixed with an original local search technique to provide better results. Numerical comparisons of solutions sets quality are made using the hypervolume metric. Based on the results for testcases that can be solved exactly, efficient implementation for PAES and NSGA2 allows execution times comparison for large instances. Results show good performances for the heuristic approaches as compared to the exact algorithm for small testcases. Approximate methods present a scalable behavior on largest problem instances. A master/slave parallelization scheme also helps to reduce execution times significantly for the modified PAES approach.
1. Introduction
Localization of facilities such as public services (schools, hospitals, etc.) is an important problem for social planners and policymakers. Most of the time, this problem is formulated as the (single objective) median problem which is a central problem in Operations Research (see, for example, [1, 2] for surveys). The median was introduced by Hakimi [3] who describes its basic properties. Its basic variant can be defined as follows: given a set of demand nodes, distance values for each pair of nodes, and a fixed number of facilities, locating each facility at one of the nodes, while minimizing the sum of distances from each node to its closest facility.
Recent developed algorithms solve the single objective version exactly for instances of thousands of nodes (e.g., 25.000 nodes in [4]). However, for a policymaker, considering additional objectives would be useful when solving median problems on real cases, leading to different variants; e.g.,(i)dispersion problem: the dispersion problem consists of spanning the facilities by maximizing the minimal distance between two of them. This objective function is suitable to locate business franchises and also when locating obnoxious facilities [5].(ii)center problem: the center problem [6] aims at minimizing the maximal distance of demand nodes from their facility, or the average distance of a fraction of those that are the farthest from their closest facility, e.g., 5% farthest of them. This problem formulation can be applied to locate emergency services such as fire stations.(iii)multimode transportation location (MTL) problem: in many real cases, transportation can be done by different means (by foot, bike, car, buses, etc.) depending on criteria like a threshold on the distance to the nearest facility. As an example, pupils are going to school by foot (category A) or by public transportation (category B), with a threshold on the distance defining the 2 categories, e.g., 2 km. The objective is then to minimize both the mean distance for those of category A and the number of people in category B.
In the following we will focus on the MTL problem with multiple objectives. The practical context is to optimize location of schools. This School Problem is a typical multimode transportation problem since, depending the distance, pupils can go to school by foot or by bus. MTL problems can be seen as multiobjective optimization problems if means of transportation have impact on each other. Optimizing the cost for one of them can degrade the other objective value. When considering a multiple objective optimization problem (MOO), a single solution optimizing all of the objectives simultaneously rarely exists. Let be an objective function mapping solutions , the search space, to the objective space , with . MOO algorithms look for solutions such that is optimized (in the sequel, we consider minimization). is a point of the objective space , with each of being one of the objective function values to be minimized. Many approaches rely on the dominance concept to choose, among a set of solutions , those ones that represent the best tradeoffs of objectives within the search space. We say that a solution dominates another one if and . It is denoted as . Namely, is as good as on all objectives and better than for at least one of them. Solutions that are not dominated by any member of are efficient solutions and constitute the Pareto set . The set of points corresponding to the efficient solutions is the Pareto front . In the sequel, we say for short that a point dominates if .
The goal of MOO algorithms is generally to determine or approximate points of and associated solutions.
Solving multiobjective median instances is of course related to median problem exact and heuristic resolution approaches but also to general approaches used for solving MOO problems, either exactly or approximately. The former are often based on Integer or Binary Programming (Multiobjective Integer Programming, MOIP), the latter on Evolutionary Multiobjective Algorithms (EMOA). Our main contribution is to develop and evaluate the two kinds of approaches, in order to be able to solve exactly medium size cases of the School Problem and approximately very large scale cases. We exploit some mathematical properties of our targeted problem in order to model it with MOIP and solve it exactly with an constraint algorithm. Large testcases are handled with 2 different general EMOA frameworks, namely, the Pareto Archived Evolution Strategy (PAES, [7]) and Nondominated Sorting Genetic Algorithm 2 (NSGA2, [8]). We have modified the former and mixed it with a local search technique: we have adapted to the multiple objective case an efficient neighborhood evaluation procedure [9] developed for the median problem. NSGA2 can also use our local search technique, as a postprocessing step. We show that, in many cases, for an equivalent computational effort, a wellknown populationbased approach such as NSGA2 is outperformed by the single individual method, PAES [7], thanks to our hybrid approach. Efficient parallelization helps for handling large testcases. As shown in Section 2, similar approaches exist for MOO median, but with different multiple objectives and local search algorithms. Furthermore, to our knowledge, no results have been presented for the parallelization of these approaches.
Section 2 introduces related work for multiple objective optimization problems embedding median like formulation. Next, in Section 3, we formalize our bicriteria multimode transportation problem, with its specific objective functions. In Section 4, we present an exact approach for solving the problem, with an constraint like technique. The problem can also be solved using popular MOO heuristic approaches like NSGA2 and PAES. We show how to adapt those frameworks to our problem solving, coupling them with an aggregation technique for performing local search. The MOO frameworks are presented in Section 5, and the local search, using limited neighborhood, is detailed in Section 5.2, along with its exploitation by the MOO methods. Evaluation and comparison of the proposed algorithms are realized with the hypervolume standard metric [10], over Beasley’s benchmark [11] in Section 6. Last, conclusion and perspectives are detailed in Section 7.
2. Related Work
A few works extend the median problem with multiple objectives. The median problem is hard [12], and MOO versions are also hard since they embed the single objective version. Thus, if some exact algorithms exist, heuristic approaches are preferred for large testcases.
The MOO median problem with an additional facility cost objective is dealt with in [13]. Each facility is weighted by a building cost, and the goal is to minimize the sum of distances to locations (median objective) and the sum of costs of the opened locations. The authors in [13] used two approaches. The first is an constraint like formulation that is a mix of twophases algorithms [14] and classical constraint approach [15]. According to the authors, it leads to a close approximation of the full Pareto front. Even if its second objective is different (facility cost instead of number of pupils by public transportation), the technique used is similar to our approach concerning the constraint problem formulation. However, since the number of nodes varies in our case for the distance count (for instance, pupils going to school by public transportation are not taken into account for the distance of demands from facilities), the method must be adapted, as it will be shown in Section 4. The second approach used in [13] helps to handle large testcases and is based on MOGA framework, with a path relinking local search procedure for mixing solutions. Problems with uniform demand at each location (unitary problems) and with up to 400 demands and 20 facilities are processed. As MOGA, NSGA2, which we are testing, also uses a populationbased approach.
In [16], the authors formulate a biobjective median problem with the following objectives : minimal distance to the closest facility (median traditional objective) and maximization of the minimal distance between facilities (dispersion problem objective). They formulate and solve it as a MOIP, with an constraint algorithm and also with a customized approach, based on threshold fixing for the minimal distance between facilities. Taking into account the threshold leads to additional constraints in a single objective subproblem. The algorithm iterates on subproblems, in order to provide the exact Pareto front of the biobjective initial problem. Once again, our second objective is different and implies an adaptation of the constraint algorithm. Furthermore, the modified MOIP developed in [16] is specific to the dispersion problem, since it iterates on a threshold value of dispersion (minimal closeness of selected locations).
Problems similar to multiobjective median have to also be examined in the context of disasters management [17], or tsunami exposure. In [18], the authors solve approximatively a 3objective optimization problem for school localization. Similarly to our heuristic method, it uses the NSGA2 framework mixed with a local search procedure, but the objective is not exactly the same as with median problems, since the number of facilities is not fixed but driven by a cost minimization objective. The second objective is related to tsunami risk exposure, and the third one takes into account both the distance to school and the number of pupils not able to go to school, as it is considered too far away (meaning, to go by foot). As we will see in the next section, this objective is in fact a mix of ours.
The median problem and its multiobjectives variants have also some military applications. Localization of repairs, supply, or security facilities can be modelized and addressed as median problems (see [19] for a survey). Common objectives are cost of localization, coverage of areas, and rapidity of deployment and also security of the localization that must be resilient. As an example, in [20], a multiobjective model, solved using a genetic annealing heuristic, takes into account cost and response time for Canadian forces prepositioning.
3. The Multimode Transportation Problem Modelling
We want to optimize the location of services, such as schools, in such a way that users go by foot to the closest service if it is possible (according to a threshold on the distance they can walk) and others take public transports. In the following, we will focus on a particular case of multimode transportation problem, the School Problem. The School Problem is stated as follows: given a set of demand nodes holding pupils and a set of candidate nodes for locating schools, a value of , and distances between each couple (demand node, candidate node) (, ), define a set of locations for setting up schools (school nodes) among the candidate nodes with the following objectives:(i)The mean distance of demand nodes to the closest school node is minimized for those demand nodes that have a distance to closest school node less than the threshold .(ii)The total number of demand nodes with a distance to the closest school greater than is minimized.
The modelling of the School Problem can be stated as follows:where(i) is the number of demand nodes;(ii) is the number of candidate nodes (number of possible locations for a school);(iii) is the number of candidate nodes to be selected as school nodes;(iv) is the threshold on walking distance;(v) is the number of pupils located at demand node ;(vi) is the distance between demand node and candidate node ;(vii) is a decision variable, indicating if the pupils at node go to school by walk or not (category A nodes).;(viii) is a decision variable, indicating if the candidate node is selected or not as a school node;(ix) is a decision variable, indicating if the demand node closest school node is or not. If , we say that (demand node) is covered by (school node) .
Equations (1) and (2) reflect the objectives stated above, taking into account the number of pupils located at each demand node. Equation (3) fixes the number of selected school nodes. Equations (4) and (5) ensure that the single selected candidate node for pupils at a demand node is effectively a school node. Equation (6) ensures that pupils at a demand node are accounted to go walking if and only if a candidate node is selected for school location in the neighborhood of the node.
This formulation induces that candidate nodes are restricted to demand nodes, but it could be extended to an arbitrary set of candidate nodes. Also notice that distance data can correspond to Euclidean distances, with known geographical locations for the different nodes, or to shortest path values, computed with an algorithm as FloydWarshall if the nodes are embedded within a routing graph, with moving cost values on the edges (e.g., time by walking, distance, and transportation cost). In the latter case, the value and meaning of the threshold are to be adapted. The use of the threshold reflects the public Swedish policy where authorities offer free transportation to some pupils based on the walking distance to school. If , we speak about the unitary version of the School Problem.
We show in the next section how to model this problem as a MOIP and how to solve it with an adhoc technique.
4. Exact Problem Formulation and Solving
It is possible to solve exactly multiple objective problems using constraint approaches. These techniques require a linear formulation of the aimed problem. Since the first objective function of our model is nonlinear (1) and nonconvex, those techniques are not directly applicable. In Section 4.1, we present the problem as a multiobjective mixed integer linear program (MOMILP). Its resolution would allow for computing the optimal value of each objective (by removing other objectives from the formulation and using a MILP solver as CPLEX, [21]). However, even computing only the mean distance (our first objective) with this model required high execution times, as detailed in Section 4.1.
Therefore, we consider a simplified MOIP in Section 4.2 where the computation of the mean distance is replaced by a computation of the cumulative distance. Using this model, we can compute the Pareto front and associated solution set using an adapted constraint [22]like approach presented in Section 4.3. We show that this method provides the exact Pareto set for the problem stated in Section 3.
4.1. MOMILP Modelling
We studied the translation of the School Problem to a MOMILP problem mainly to compute directly the minimum mean distance (1) using a MILP solver. Since this objective function is continuous, its linearization requires using continuous variables. Indeed, starting from the modelling of the School Problem presented in Section 3, we can linearize (1) using a new set of (continuous) variables , which replaces the set of (Boolean) variables . The semantics of variables is the following: if demand node is covered by school node and is in category (i.e., ) and, in this case, .
Using this semantics, objective (1) becomes linear, as shown in . Notice that, by removing the variables , we “forget”, for each demand node in category (i.e., with ), which school node covers it. It is possible to keep this information, but it is useless for the resolution of the problem (as long as there is at least one school).
To satisfy the semantics of the variables and specifically specify , we also add a continuous variable for each demand node and a continuous variable . The resulting model is as follows:
In this model, , and , ensure that : from and (and since are binaries), we have . By , we get , which gives the lower bound of .
From this result, states that when demand node is in category (i.e., ), (the minimization of ensures that only one is equal to and the other are zero). Equation implements Constraint (4), whereas the implementation of Constraint (6) is done by and (note that the former is needed for the first objective while the latter is needed for the second one).
We have tried to exploit this formulation for computing the minimal mean distance, with the CPLEX MILP solver [21], using the single objective function . In practice this approach did not work, except on our smallest example (100 nodes and 5 schools). On the other testcases, CPLEX returned possibly nonoptimal solutions after reaching its memory limit. We think that this failure stems from two issues. The first issue is the variability of (and, as a result, the values of all nonzero ) which makes node selection a hard problem. The second issue is that when solving relaxed problems, with and partly continuous, enables putting even when is close to 1, since is very small. As a result, the algorithm needs to force almost all (and ) to integral values before getting a meaningful minimal bound. Since both issues are related to the computation of a mean distance, we considered a MOIP using the cumulative distance instead.
4.2. MOIP Modelling
The MOIP formulation of the School Problem with the cumulative distance can be seen as a simplification of the MOMILP formulation of the previous section. In fact, we just need to state that . As a result, the variables are useless (, , and disappear). From this result, we propose the MOIP formulation for the School Problem, as follows:
Compared with the previous model, we replace by (12), (13), and (14), which ensure that demand node cannot be in if there is no school node at distance .
Notice that (13) and (15) both use , but with different meanings: the former ensures that a pupil cannot walk to a distant node (to minimize the second objective), while the latter prevents an artificial minimization of the cumulative distance by covering a demand by a distant node where a closer one is available.
As we have seen before, our School Problem is formalized only approximately by , since (9) uses the cumulative distance instead of the mean distance for pupils in category A. While the minimization of the cumulative distance can be compared to the median problem (indeed, we can get the median problem by adding for all and removing (13)), its exact resolution using MOIP solvers can be more difficult, as the relaxed problem (with continuous variables) has more solutions. To explain this issue, let us consider only two demand nodes (which are also candidate nodes) 1 and 2 such that and one location (). With the relaxed median problem (), the solver gives directly the optimal cumulative distance , whatever the values of and . Using the relaxation of , the optimal solution sets all variable to except and which are set to 0 and the cumulative distance is 0.
Using this observation and the fact that the cumulative distance is still not the mean distance lead us to consider an adapted constraint method where the first objective is (10).
4.3. Exact Algorithm
constraint [22] methods proceed as follows: a series of single objective (i.e., Integer Program, IP) problems are solved, transforming all but one of the objectives of the MOIP considered into IP constraints. The constraint set is updated at each iteration to enforce the exploration of the whole objective space. In their paper [23], Ozlen and Azizoğlu introduce a recursive algorithm to generate a Pareto set for a MOIP problem. They use the set of already solved subproblems and their solutions to avoid solving a large number of IPs. In [14], a two phases algorithm is applied for the biobjective assignment problem. In the 2 objectives case, it first determines the extreme points in the objective space by discarding one objective at a time and solving the resulting single objective problem. Then, in a second phase, it partitions the objective space according to the range of values found at the first step and explores it by slices. It also combines the approach with heuristics specific to the assignment problem for enhancing the execution times.
An adapted constraint algorithm which computes all nondominated points sorted by values (i.e., second objective) is presented in Algorithm 1 for solving . This approach uses the fact that given a nondominated point , other nondominated point such that must satisfy both and . Since the minimal cumulative distance is hard to compute (and not really useful as it may not be the minimal mean distance), we do not compute the extreme point associated with it.

In , we consider the minimization of the cumulative distance instead of the mean distance. Hence both objective function parameters are of integer values (once all distances are converted into integers), which ensures that our constraint method cannot miss any efficient solution for the cumulative distance problem. Furthermore, one can check that, for any feasible solution of value , the mean distance for people in category is where . Especially if for all (Unitary Problem). Hence, any efficient solution for the mean distance problem is also an efficient solution of the cumulative distance problem.
Theorem 1. Let be a feasible solution for (associated with point in objective space), such that it is efficient for mean distance optimization, i.e., for each feasible solution (point ):Then is efficient for cumulative distance optimization.
Proof. If , the result is straightforward. Otherwise, and , which implies .
Thus all efficient solutions for the School Problem are found by solving . But the converse is not true, and nondominated points can also appear in the Pareto front for that do not belong to the Pareto front for the School Problem. However, since the solutions generated are sorted by the value of , it is easy to make the algorithm generate only nondominated points for the School Problem: for each generated nondominated point , we add the following constraint to the problem:
Constraint (21) (which has only integer coefficients) ensures that any subsequent solution has a mean distance strictly less than the previous one. Including this constraint, Algorithm 1 computes the Pareto front for the School Problem.
5. Heuristic Problem Solving
MOO problems search spaces are often intractably large [15]. Many heuristics have been developed for searching over huge spaces, particularly using Evolutionary Multiobjective Algorithms (EMOA). Genetic algorithms have been shown to be interesting for solving very large instances of the median problem in [24]. We investigate the extension of this work to the MOO case, focusing on chromosomebased EMOA. We investigate the use of two algorithms, NSGA2 and PAES.
NSGA2 is a reference algorithm in the EMOA field, and it compares positively to PAES for standard benchmark problems [8]. But as stated below, NSGA2 (and many other EMOAs) iterates onto a population of individuals (each of them represents a solution), producing eventually a large offspring at each generation. Even if this is an advantage when exploring the search space for building an interesting front, it can be very costly when employing such an algorithm for solving a problem with a costly evaluation function. The evaluation function can be intrinsically complicated (e.g., physics simulation) or time costly because it runs a local search within the neighborhood of the solution to be evaluated. We have already applied PAES successfully for the former reason in the field of RealTime Scheduling [25]. We have also used a local search technique for the single objective median problem in [24]. Thus our goal is to compare PAES coupled with a local search technique to the reference algorithm NSGA2 for the MOO School Problem.
We first present and compare shortly NSGA2 and PAES and then we detail our declination of both frameworks with a local search technique for solving the School Problem.
5.1. PAES and NSGA2
In many evolutionary algorithms, a solution is represented by a set of parameters called chromosome (or genotype). The encoding of solutions (i.e., the data structure of a chromosome) is specific for each problem. Several researchers have proposed genetic algorithms (GA) for solving the median problem [26, 27]. Most of them use a classical string representation; i.e., each chromosome is represented by a single array of integers of length embedding the index of the selected candidate nodes. As stated in [24], we will use the same encoding, i.e., the chromosome represents the list of school nodes in our case. We add the constraint that in all chromosomes no school is duplicated. The initial population (or first solution) of our algorithms is generated by picking up distinct random candidate nodes, leading to feasible solutions. All the candidate nodes have the same probability to be chosen.
Both EMOA’s goals are to produce a front that preserves diversity of solution objective values, with points that are numerous and spread well along and the accuracy of this front, with points that are close to . Within a population set of solutions, both PAES and NSGA2 use a crowding metric for handling diversity: PAES by defining an adaptive grid over the objective space and counting points within each portion of the grid and NSGA2 by measuring the distance between points and their closest neighbor in each direction for each objective. To ensure accuracy, both algorithms embed an elitism mechanism: the number of solutions that are kept across generations (elitism) is given either by the population size (NSGA2) or by an external archive (PAES). This size also controls the number aspect of diversity. In order to increase accuracy of solutions and convergence to , PAES uses dominance locally, comparing current solution to a single offspring at each generation, while NSGA2 sorts its population into dominance classes. PAES eventually updates its archive with the offspring, while NSGA2 renews its population by generating a series of offsprings at each generation, with a dominance class based selection. With both algorithms, when comparing solutions, if they are nondominated against each other, the crowding metric is used to tie them. Concerning the operators used for generating offspring, we have implemented a specific mutation operator that preserves feasibility within PAES and NSGA2 applies standard binary mutation and crossover onto a binary representation of chromosomes it generates internally. We also provide a penalty metric that helps NSGA2 discarding unfeasible solutions (penalizing duplicated schools).
We will see in the two next sections how to adapt those strategies for mixing the various approaches with a fast local search technique inspired from median single objective literature.
5.2. Local Search Technique for MOO
Many heuristic methods have been developed for the (single objective) median problem (see surveys [1, 2]). Variable neighborhood search [9] is very popular, coupled with fast evaluation techniques [1]. Those techniques, based on precomputing of the 2 closest schools for each demand node, allow updating the solution cost faster when modifying only partially the solution itself during the local search process: for each demand node, for a given set of schools (solution), the distance to closest school (, the one that covers the demand node) and the second closest school () and its distance to the node are stored. We use a neighborhood operator called hypermutation, inspired from [24]. The neighborhood size is controlled by the number of modified school nodes (seeds): the neighborhood of a solution, for a fixed number of seeds , is defined by replacing one by one each seed node by all of the possible other candidate nodes; i.e., those that are not yet selected as school nodes within the solution. The solution cost can be updated faster, using the precomputed tables mentioned above: when a school is replaced by , demand nodes covered by it () are now covered by their second closest school (), except if is closer to them than . Objective functions are updated according to the case. We fix the number of seeds to for a School Problem, to obtain a reasonable neighborhood size. This size is , with the number of seeds, the number of candidate nodes, and the number of schools.
The neighborhood and associated evaluation techniques are devoted to the single objective median problem. They can be kept when running the MOO case. However, evaluating and comparing to current front multiple neighbors resulting from multiple function evaluations could be very time consuming: evaluation itself can be costly if applied to a large number of neighbors, and these numerous solutions must by compared against existing population or archive, with a dominance checking procedure. Thus, we propose an approach where the neighborhood exploration is run with a single objective metric and only the most interesting solutions found during the local search are compared to the current solutions kept in population or archive. This algorithm is outlined in Algorithm 2. We transform, for the local exploration, the 2 objectives metric of the School Problem into a single one by a classical aggregation technique: weights are given to objectives and is computed. and are normalized values (according to the objective vector of the neighborhood search starting point). is evaluated for different values of in order to explore the objective space into different directions (e.g., ). Since is restricted to a few values, only a subset of the supported solutions can be found and of course unsupported ones cannot be detected by the local search [28]. For each direction searched, only the best solution (i.e., with the minimal value for the associated weight) into the neighborhood is provided as a result at the end of the local exploration process. This approach can be related to the decomposition techniques used in algorithms like MOEA/D [29].

5.2.1. Cost of Evaluation
The grain of the search is tuned by the number of seeds , leading to a neighborhood of size , and also by the number of different values for . Let by the cost of a (standard) evaluation of a solution. It mainly consists of finding the closest school for each children and is in . Evaluating a series of neighbors for a neighborhood of () costs , also stated as . During the (fast) evaluation of a solution as a neighbor of another one, only a single school is modified and only the distances to school for its covered nodes are to be updated, in time. Thus evaluating a whole neighborhood takes , with roughly, a factor of gain as compared to .
5.3. Mixing the Local Search Technique with EMOA Algorithms
This local exploration procedure can be mixed with the different EMOAs into 2 different ways:(i)Applying the local search for each individual: the offspring of an individual is generated via the local search procedure, with directions of search.(ii)Using the local search as a refinement step: applying local search to the members of the population or archive obtained as a postprocessing step.
NSGA2 processes natively multiple offsprings, but the original PAES manipulates a single current solution and a single offspring (it is a procedure) that replaces (or not) the current solution at next generation, depending on dominance checking. In the first way of mixing, a algorithm such as PAES must be adapted. Knowles proposes such a version in [7]. We use instead the same selection mechanism as in [30], developed for costly functions coupling with PAES and well adapted in the context of an hybrid method combining EMOA with (time costly) local search. The selection procedure is as follows: the offsprings resulting from the local search are compared to the PAES archive using PAES policy. Then, the new current solution is chosen considering the whole archive as follows, based on a randomly selected optimization criteria: a random integer in the interval is generated. If it corresponds to an objective function (i.e., ), these criteria are chosen; otherwise () the crowding criteria is selected. The next current individual is chosen randomly within the 10% best ones within the archive for the elected criteria. Concerning crowding, the best individuals are those in the less crowded areas of the grid.
The goal of the approach is to solve large size problems with a high execution time cost induced by evaluation and proportionally a small amount of time devoted to EMOA internal computations. So we consider here that the cost of choosing next current individual by processing the whole archive is not a drawback for execution times according to the benefit in terms of quality of search that we expect.
6. Algorithms Comparison
The multimode transportation problem or School Problem to be solved is specific, and to our knowledge, no comparable algorithm exists for solving it (see Section 2). Thus, this evaluation section mainly compares the heuristic approaches that we have developed to our customized constraint method. The following algorithms are compared:(i)Exact: the constraint method presented in Section 4, which provides the exact Pareto Front ,(ii)EMOA approaches presented in Section 5, using NSGA2 and PAES algorithms. For both EMOAs, the standard and the hybrid versions tested: local search realized at each generation with PAES and as a postprocessing step for NSGA2.
We first describe the benchmark set used for comparing algorithms results, detail algorithms setup, and then present results for both quality and execution time aspects.
(a) Data Set. The data set used for the evaluation is the Beasley benchmark set, devoted to the median problem [11] and widely used for testing median solvers [26, 31, 32]. It is a set of 40 testcases. Candidate and demand nodes are identical, with a unitary demand at each node. (the number of nodes) varies from 100 to 900 and value from 5 to 200. Testcases are ordered by increasing value and, for a given , by increasing value as shown in Table 1.

We compute the threshold in such a way that 15% of the distances between nodes are lower than .
The largest graph of this benchmark has 900 candidate nodes and is not particularly large according to the size of problems current median solvers can manage (e.g., instances of up to 24,978 nodes solved exactly in [4]). However, multiple objective optimization could lead to much more costly evaluation of the solution space for the same instances as compared to singleobjective formulations, and Beasley benchmark contains large enough instances for our purpose, as execution times will show. It is also comparable to other biobjective median problem instances used with others approaches (e.g., 402 nodes for the largest instance in [13]).
(b) Algorithms Setup. The exact method has been implemented in C. It is a modified version of the AIRA software, a general purpose MOIP solver [33], and it uses CPLEX solver 12.3 as IP solver software [21]. As heuristic, we use the version of NSGA2 provided by Deb at [34], with the following parameters: a population size of 100, and 500 generations. We also use our own implementation of PAES, based on the C version provided by J. Knowles at [35]. The depth parameter of the algorithm is used to define the grain of the grid used for measuring the crowding of the objective space by the members of the PAES archive. It defines the number of recursive subdivisions of the range of values for each objective. It is set to 4, according to PAES recommendations for some biobjective problems. The archive size for PAES is of 100 and the number of generations of 1000. The parameters for the local search are (number of directions for search), corresponding to and for the number of seeds, since is always greater than for the Beasley benchmark testcases (see Table 1). Those parameters have been set in order to equilibrate the computational effort of the different algorithms, as execution times comparison will show. They lead to 1000 (resp. 50.000) standard evaluations and 25.000 (resp. 2.500) fast evaluations for hybrid PAES (resp. hybrid NSGA2) (see Section 5.2). All of the tests are performed onto a 48 processors (Xeon E5 at 2.2GHz) SMP machine, with 132GB of memory, and running Linux CentOS.
(c) Quality Evaluation Procedure. The sets of solutions provided by the algorithms are to be compared qualitatively. Many metrics exist [28] that must take into account both the quality of solutions obtained (accuracy) and their number and location along the Pareto front range (diversity). The hypervolume (HV) [10] is often used to compare 2 sets of solutions and . It computes the area defined by each set; according to a reference point dominated by all of the points in the sets, the highest value is the best one. Figure 1 illustrates the comparison by HV for 2 fronts, with two objective functions and to be minimized.
For each testcase, we run each algorithm in competition 30 times. We collect all of the resulting individuals and compute the associated nondominated front. This front’s bounds on objectives are then computed and used for normalizing the objective vectors associated with the individuals into the range . The point is thus dominated by all of the vectors and can be used for computing hypervolumes (see values in Figure 1), leading to a maximal possible HV of , if a single objective vector dominates all of the others. We provide results for each algorithm and also make a statistical analysis with the KruskalWallis nonparametric test [37], for comparing the quality value series obtained by running each algorithm 30 times for a given testcase: if and only if a first test for significance of any differences between the samples is passed (H0 hypothesis), at a given confidence value (we use the standard 0.05 value), then the output will be the onetailored p values for rejecting a null hypothesis of no significant difference between two samples, for each pairwise combination. If the p value for a testcase is less than 0.05, we considerate that one algorithm beats the other for this testcase with a sufficient confidence. On the contrary, if series are not different enough, or the p value obtained for characterizing the differences is over 0.05, the testcase is not taken into account for comparing the couple of algorithms.
6.1. Hypervolumes Comparison
We compare the hypervolumes obtained for 5 algorithms:(i)the Exact algorithm described in Section 4. The 10 first testcases of the Beasley benchmark can be solved in less than 1 hour with this method, providing a reference front for comparison with heuristics for this subset of the Beasley benchmark.(ii)the PAES algorithm, with our selection method, but without local search.(iii)the hybrid PAES algorithm, embedding local search at each generation.(iv)the NSGA2 algorithm,(v)and the hybrid NSGA2 algorithm, using the local search as a postoptimization method
All the heuristics are applied onto the whole benchmark, with parameters as previously described. We have removed from the resulting mean values 5 testcases (#20, #24, #25, #30, and #34): for those ones, NSGA2 (and thus also hybrid NSGA2) does not provide any feasible solution for at least 5 runs (all runs for #30). This is due to the fact that unfeasibility corresponds to duplicated schools in solutions, and this can arise with a higher probability when increases. For the discarded testcases, is between 100 and 200 (see Table 1). Statistics are thus calculated over the 35 remaining testcases.
Figure 2 shows the average values of hypervolumes for the different algorithms in competition. On average over the runs for the 40 testcases, hybrid PAES provides the best results for all testcases, when considering only EMOAs. It outperforms, respectively, PAES of 24.9%, NSGA2 of 58.7%, and hybrid NSGA2 of 48.0% on average. For the 10 testcases for which exact solution is known, it provides a front with an HV that is, on average, at 7% of the optimal one. Clearly, the generic NSGA2 framework do not allow to obtain good results, as compared to the specialization of PAES for the School Problem: within NSGA2 version, the only specific component is the objective function, with nonspecialized mutation and crossover operators that can lead to unfeasible solutions. This is managed with NSGA2 by constraints penalties, but, unfortunately, this mechanism does not allow an efficient search of the feasible solution space. As a symptom, the size of the fronts provided by the NSGA2 algorithm is on average 68% less that those found by hybrid PAES (63% less for hybrid NSGA2). It is 49% when comparing PAES and hybrid PAES, showing that local search effectively helps in discovering nondominated solutions. The result is that local search allows improving results: hybrid PAES outperforms PAES by 24.9% and the local search by hybrid NSGA2 allows improving NSGA2 results in terms of HV by 32.7% on average. The latter is very significant, for a single local search phase, but the improvement is realized on the (relatively to hybrid PAES) poor results obtained with NSGA2, which leads room for optimization.
Table 2 assesses the confidence that can be given when comparing those mean results, by applying the KruskalWallis nonparametric test to the hypervolume series. Even if hybrid PAES always provides mean HV values better than other EMOAs in competition, this is not completely assessed by the statistical test for all of them; e.g., hybrid PAES beats hybrid NSGA2 36 times only with the defined confidence of 0.05. This can be explained by the closeness of the results in some cases: for example, for testcase #6, mean HV is of 0.4061 for hybrid PAES and 0.3937 for hybrid NSGA2, with respective standard deviation of 0.002 and 0.030, meaning that result values of the two algorithms overlap.

6.2. Execution Times
Average execution times of the different heuristic approaches are depicted in Figure 3 for testcases of Figure 2.
PAES algorithm runs in less than a second for all of the testcases, thanks to the single evaluation at each generation. Because of the local search executed only as a postprocessing step, hybrid NSGA2 execution time is very close to NSGA2’s one, with a local search realized, on average, in 1.62 seconds for the 100 individuals of the population. One can see that the same order of computational effort has been used for the different algorithms, by tuning population size and number of generations ( and for NSGA2 versions, for PAES and its hybrid version). Those parameters allow for exhibition experimentally that execution times of tested algorithms depend on a combination of and (as shown in Section 5.2). Reminding the characteristics of the problems in Table 1, NSGA2 (and the hybrid version) beats hybrid PAES for values of and if the number of nodes , and hybrid PAES is faster in the other cases. Furthermore, NSGA2 execution times are positively correlated to the chromosome size for genetic operations, i.e., to the value of (e.g., for instances 1 to 5, with from 5 to 33). The decreasing steps in curve for the hybrid PAES correspond also to the increasing values of , for a fixed number of nodes (e.g., for instances 26 to 30), but with the inverse effect as the one depicted for NSGA2. Concerning hybrid PAES, for the local search, the number of seeds is fixed (see Section 5.2), the number of alternatives for each seed is of (), leading to a size of for the neighborhood. Thus, the more the value, the smaller the neighborhood.
(a) Parallelization. Hybrid PAES is promising in terms of quality as compared to (hybrid) NSGA2 but time consuming as compared to PAES. We have implemented a masterslave scheme for realizing the evaluation and local search in parallel for different individuals, inspired from [30]. Parallel hybrid PAES implementation (in C) uses multiple threads (one per slave plus one for the master process). The OS scheduler allocates each of them to a core if the number of cores is sufficient. Figure 4 shows the speedups obtained when running the algorithm with more than one slave (and thus more than 2 cores) for the evaluation and local search of individuals. The execution times are obtained using a 48 cores SMP system, ensuring that each slave (and the master) is run on a separate core. In order to increase the computational effort, the number of generations is set to 5000 instead of 1000 in previous tests.
Clearly, adding computing resources helps in speeding up the computations and especially when the computational effort is high (last testcases). The average speedup is 1.89 (resp. 3.76, 7.48, 14.07, and 18.10) for 2 (resp. 4, 8, 16, and 32) slaves. Speedups look linear according to the number of slaves for up to 8 slaves. Classically, the cost of the parallelization is due to the bottleneck in communications induced by the masterslave scheme and also to parallelization itself: for example, for testcase #35, the sequential version of hybrid PAES runs in 29 seconds for 1000 runs (see Figure 3) and the parallel version, with a single slave, runs in 273 seconds for 5000 generations (instead of expected 29 seconds 5 = 145 seconds). The communication bottleneck degrades the speedup for 16 slaves, but execution times are still improved in all cases as compared to 8 slaves version, except for testcase #5 (average speedup of 14.07). But for 32 slaves, the improvement of speedup is very limited (18.10 on average). As we find again the same steps shape on curves as those observed in Figure 3 for sequential times, with roughly higher speedups for the most costly testcases, it is possible that the parallelism grain is not sufficient to ensure good acceleration with 32 slaves. Notice that this can be tuned by both (number of directions for local search) and (number of seeds) parameters of the hybrid PAES (see Section 5.2).
Parallelization helps in reducing the drawback of large execution times for hybrid PAES. For example, the largest execution time (testcase #35) is 273 seconds for the algorithm with a single slave and 5000 generations, and it is decreased to 13 seconds with 32 slaves. Since speedups increase with both the number of slaves and the size of the problem handled, the parallel approach for hybrid PAES looks promising for the processing of real very large data.
7. Conclusion and Future Work
We present in this paper a MOO modelling of a facilities localization problem, applied to a multimode transportation problem, the School Problem. The goal is to optimize the transportation mode for pupils, with two possible modes, namely, by foot and by public transport, with constraints on the walking solution. A modelling for the School Problem is proposed, and an exact resolution method, based onto an Integer Programming formulation, is defined. The IPbased approach solves in a few minutes to one hour each of the 10 smallest instances of the Beasley public benchmark. A heuristic method called hybrid PAES mixes median neighborhood search with PAES EMOA. Our approach is able to solve the same 10 instances in less than 10 seconds each, with an average degradation of quality (measured with the hypervolume metric) of 7%. Hybrid PAES also provides solutions for large instances for which the Pareto front is unknown. For those cases, it is competitive with a standard approach as NSGA2, with results that outperform this reference method by 58% and 49% for hybrid NSGA2. The execution times are also improved as compared to NSGA2 when the problem complexity is significant enough ( and ). A parallel version of hybrid PAES is also proposed, using a masterslave scheme. The speedup of the algorithm is linear for up to 8 slaves, close to 14 for 16 slaves, but degrades with more slaves (speedup is only 18.2 for 32 slaves). We think that this limitation on the parallelization degree is due to the grain of parallelism induced by the data and the algorithm’s setup for the experiments and that it will not hold for larger instances. We plan to apply it to real data, as realized in [24] for single objective median problems, dealing with countryscale instances, with nonuniform population distribution. This is important because real problems may exhibit particularities in their data. Another direction of research is to add a third objective function, related to facility implantation costs, with a relaxation of the number of facilities ( becomes a bound instead of a fixed value). This economic cost is useful, for example, in disasters management.
Conflicts of Interest
The authors declare that there is are conflicts of interest regarding the publication of this article.
References
 N. Mladenovi, J. Brimberg, P. Hansen, and J. A. MorenoPrez, “The pmedian problem: A survey of metaheuristic approaches,” European Journal of Operational Research, vol. 179, no. 3, pp. 927–939, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 J. Reese, “Solution methods for the pmedian problem: An annotated bibliography,” Networks, vol. 48, no. 3, pp. 125–142, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 S. L. Hakimi, “Optimum locations of switching centers and the absolute centers and medians of a graph,” Operations Research, vol. 12, no. 3, pp. 450–459, 1964. View at: Publisher Site  Google Scholar
 S. Garca, M. Labbé, and A. Marn, “Solving large pmedian problems with a radius formulation,” INFORMS Journal on Computing, vol. 23, no. 4, pp. 546–556, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 E. Erkut, “The discrete pdispersion problem,” European Journal of Operational Research, vol. 46, no. 1, pp. 48–60, 1990. View at: Publisher Site  Google Scholar  MathSciNet
 N. Mladenović, M. Labbé, and P. Hansen, “Solving the pcenter problem with tabu search and variable neighborhood search,” Networks, vol. 42, no. 1, pp. 48–64, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 J. D. Knowles and D. W. Corne, “Approximating the nondominated front using the pareto archived evolution strategy,” Evolutionary Computation, vol. 8, no. 2, pp. 149–172, 2000. View at: Publisher Site  Google Scholar
 K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan, “A fast elitist nondominated sorting genetic algorithm for multiobjective optimization: NSGAII,” in Parallel Problem Solving from Nature PPSN VI, vol. 1917 of Lecture Notes in Computer Science, pp. 849–858, Springer, Berlin, Germany, 2000. View at: Publisher Site  Google Scholar
 P. Hansen and N. Mladenovi, “Variable neighborhood search: principles and applications,” European Journal of Operational Research, vol. 130, no. 3, pp. 449–467, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 C. M. Fonseca, L. Paquete, and M. LópezIbáñez, “An improved dimensionsweep algorithm for the hypervolume indicator,” in Proceedings of the IEEE Congress on Evolutionary Computation (CEC '06), pp. 1157–1163, IEEE Press, July 2006. View at: Google Scholar
 J. E. Beasley, “ORLibrary: distributing test problems by electronic mail,” Journal of the Operational Research Society, vol. 41, no. 11, pp. 1069–1072, 1990. View at: Publisher Site  Google Scholar
 O. Kariv and S. L. Hakimi, “An algorithmic approach to network location problems II: the pmedians,” SIAM Journal on Applied Mathematics, vol. 37, no. 3, pp. 539–560, 1979. View at: Publisher Site  Google Scholar
 J. E. C. Arroyo, P. M. Dos Santos, M. S. Soares, and A. G. Santos, “A multiobjective genetic algorithm with path relinking for the pmedian problem,” in Advances in Artificial Intelligence IBERAMIA 2010, A. KuriMorales and G. Simari, Eds., vol. 6433 of Lecture Notes in Computer Science, pp. 70–79, 2010. View at: Publisher Site  Google Scholar
 A. Przybylski, X. Gandibleux, and M. Ehrgott, “Two phase algorithms for the biobjective assignment problem,” European Journal of Operational Research, vol. 185, no. 2, pp. 509–533, 2008. View at: Publisher Site  Google Scholar  MathSciNet
 C. A. C. Coello, G. B. Lamont, and D. A. van Veldhuizen, Evolutionary Algorithms for Solving MultiObjective Problems, Springer, New York, NY, USA, 2007. View at: MathSciNet
 F. Sayyady, G. K. Tutunchi, and Y. Fathi, “Pmedian and pdispersion problems: A bicriteria analysis,” Computers & Operations Research, vol. 61, pp. 46–55, 2015. View at: Publisher Site  Google Scholar
 R. Abounacer, M. Rekik, and J. Renaud, “An exact solution approach for multiobjective locationtransportation problem for disaster response,” Computers & Operations Research, vol. 41, no. 1, pp. 83–93, 2014. View at: Publisher Site  Google Scholar
 K. F. Doerner, W. J. Gutjahr, and P. C. Nolz, “Multicriteria location planning for public facilities in tsunamiprone coastal areas,” OR Spectrum, vol. 31, no. 3, pp. 651–678, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 J. E. Bell and S. E. Griffis, “Military applications of location analysis,” Applications of Location Analysis, vol. 232, pp. 403–433, 2015. View at: Publisher Site  Google Scholar
 A. Ghanmi and R. H. A. D. Shaw, “Modelling and analysis of Canadian Forces strategic lift and prepositioning options,” Journal of the Operational Research Society, vol. 59, no. 12, pp. 1591–1602, 2008. View at: Publisher Site  Google Scholar
 “CPLEX online reference manual,” https://www.ibm.com/support/knowledgecenter/SSSA5P_12.7.1/ilog.odms.studio.help/Optimization_Studio/topics/COS_home.html. View at: Google Scholar
 J.F. Bérubé, M. Gendreau, and J.Y. Potvin, “An exact ϵconstraint method for biobjective combinatorial optimization problems: application to the traveling salesman problem with profits,” European Journal of Operational Research, vol. 194, no. 1, pp. 39–50, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 M. Ozlen, B. A. Burton, and C. A. MacRae, “Multiobjective integer programming: an improved recursive algorithm,” Journal of Optimization Theory and Applications, vol. 160, no. 2, pp. 470–482, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 P. Rebreyend, L. Lemarchand, and R. Euler, “A computational comparison of different algorithms for very large pmedian problems,” in Proceedings of the 15th European Conference on Evolutionary Computation in Combinatorial Optimization, vol. 9026 of Evolutionary Computation in Combinatorial Optimization, pp. 13–24, Springer International Publishing, Copenhagen, Denmark, April 2015. View at: Google Scholar
 R. Bouaziz, L. Lemarchand, F. Singhoff, B. Zalila, and M. Jmaiel, “Architecture exploration of realtime systems based on multiobjective optimization,” in Proceedings of the 20th International Conference on Engineering of Complex Computer Systems, ICECCS 2015, pp. 1–10, Gold Coast, QLD, Australia, December 2015. View at: Publisher Site  Google Scholar
 E. Correa, M. Steiner, A. A. Freitas, and C. Carieri, “A genetic algorithm for the pmedian problem,” in Proceedings of the 2001 Genetic and Evolutionary Computation Conference (GECCO2001), pp. 1268–1275, 2001. View at: Google Scholar
 O. Alp, E. Erkut, and Z. Drezner, “An efficient genetic algorithm for the pmedian problem,” Annals of Operations Research, vol. 122, pp. 21–42, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 C. A. C. Coello, G. B. Lamont, and D. A. van Veldhuizen, Evolutionary Algorithms for Solving MultiObjective Problems, SpringerVerlag, New York, NY, USA, 2nd edition, 2007. View at: MathSciNet
 Q. Zhang and H. Li, “MOEA/D: a multiobjective evolutionary algorithm based on decomposition,” IEEE Transactions on Evolutionary Computation, vol. 11, no. 6, pp. 712–731, 2007. View at: Publisher Site  Google Scholar
 R. Bouaziz, L. Lemarchand, F. Singhoff, B. Zalila, and M. Jmaiel, “Efficient parallel multiobjective optimization for realtime systems software design exploration,” in Proceedings of the 2016 27th International Symposium on Rapid System Prototyping: Shortening the Path from Specification to Prototype, RSP 2016, pp. 58–64, USA, October 2016. View at: Publisher Site  Google Scholar
 J.C. Gay, Résolution du problãme du pmédian, application ã la restructuration de bases de données semistructurées [Ph.D. thesis], Université BlaisePascal, ClermontII, France, 2011.
 A. Alkhedhairi, “Simulated annealing metaheuristic for solving pmedian problem,” International Journal of Contemporary Mathematical Sciences, vol. 3, no. 2528, pp. 1357–1365, 2008. View at: Google Scholar  MathSciNet
 “AIRA software,” Available: https://bitbucket.org/melihozlen/moip_aira. View at: Google Scholar
 “NSGA2 software,” Available: http://www.egr.msu.edu/kdeb/codes.shtml. View at: Google Scholar
 “PAES software,” Available: https://www.cs.bham.ac.uk/~jdk/mult. View at: Google Scholar
 R. Bouaziz, L. Lemarchand, F. Singhoff, B. Zalila, and M. Jmaiel, “Multiobjective design exploration approach for ravenscar realtime systems,” RealTime Systems, vol. 54, no. 2, pp. 424–483, 2018. View at: Publisher Site  Google Scholar
 W. Conover, Practical Nonparametric Statistics, Wiley, 3rd edition, 1999.
Copyright
Copyright © 2018 Laurent Lemarchand et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.