1. Introduction
The
p-median problem is one of the most well-studied facility location problems in the field of management and operation research. In the classical
p-median problem, a certain number of
p facilities are located among
n candidate locations and
m demand points are allocated to the
p opened facilities, while the objective is to minimize the sum of the transportation cost/distance to serve all of the demand points. The
p-median problem and its extensions have been widely used in many practical situations including the location of plants, warehouses, distribution centers, hubs, and public service facilities [
1]. Previous study has proven that the
p-median problem is Non-deterministic Polynomial (NP) hard [
2], and that its optimal solution is unlikely to be obtained through an exact approach with polynomial complexity. Thus, various approximate heuristics have been developed to find optimal or near optimal solutions.
The genetic algorithm (GA) is a kind of metaheuristic inspired by the biological evolution process, which is one of the most efficient metaheuristics for solving the
p-median problem. Hosage and Goodchild [
3] first reported the application of the GA for the
p-median problem, and showed the advantage of the GA in the generality for approaching difficult optimization problems. Alp et al. [
4] improved the efficiency of the GA for solving the
p-median problem by integrating the greedy heuristic in the evolution process, which showed better performance than the simulated annealing algorithm proposed by Chiyoshi and Galvão [
5]. Correa et al. [
6] employed a heuristic “hypermutation” operator in the GA to solve the capacitated
p-median problem, which obtained better results from a comparison with the Tabu search algorithm. Alba and Dominguez [
7] compared a class of GA based heuristics with several other metaheuristics and reported their relative performances. Their experimental results indicated that the GA based algorithm had advantages of wide applicability, low implementation effort, and high accuracy. Kratica et al. [
8] proposed two GAs with new coding schemes and genetic operators for solving the single allocation
p-hub median problem. Li et al. [
9] investigated the impacts of initial population generating strategies on the performance of the GA, and proposed an approach to enhance the quality of the initial population. Experimental results showed that the computational time could be reduced by improving the initial population. Herda [
10] introduced the parallel strategy in the GA to improve the efficiency in solving the capacitated
p-median problem, and designed two parallel GAs that could run on computing clusters with high performance. Memari et al. [
11] proposed a modified firefly algorithm, which was used to solve a supply chain network problem. GA has also been used to efficiently solve many other problems in the industrial field such as material handling control [
12], the transportation–production problem [
13], and the works transport organization and control problem [
14]. More metaheuristics in the facility location field can be found in Mladenović et al. [
1].
The
p-median facility location problem is usually a long-term strategic decision problem, and most of the above studies assume that the facilities’ operational environments is deterministic. However, in practical situations, the fierce competition, changes in customer taste, community growth, and economic vitality in the market area all make the demand quite uncertain in the planning period and difficult to estimate. In this situation, the demand is assumed to be an uncertain parameter, which can be modeled as a set of alternative future scenarios. Then, robust optimization approaches can be utilized to minimize the maximum regret of the overall transportation cost [
15]. Weaver and Church [
16] studied the
p-median problem in a stochastic network, where the demand in each node was modeled as a set of discrete scenarios. Chen et al. [
17] investigated the uncertain
p-median problem, and developed an
α-reliable mean-excess model to minimize the expected regret. Snyder and Daskin [
18] introduced the
β-robust approach for the
p-median problem with uncertainties, which can be used to obtain a trade-off between the expected transportation cost/distance and the robustness. The
β-robust approach for facility location integrates the advantages of both the stochastic and robust optimization methods. Then, Lim and Sonmez [
19] further studied the
β-robust facility relocation problem under budget constraint. Sakalli and Atabas [
20] employed both ant colony optimization and the GA to solve a tactical production–distribution planning problem in an uncertain environment. In most of the studies on robust facility location, in order to estimate the regrets under different scenarios, the optimal solutions for all scenarios should be first obtained. The instance under each scenario usually corresponds to a kind of classical facility location problem, which is NP-hard. The solving of each scenario problem has to employ approximate heuristics such as the GA and Lagrangian heuristic [
21], which is time-consuming. When the number of scenarios is lower, all of the instances can be solved through iteratively running an approximate heuristic. For example, nine scenarios were included in the experiments of both [
18,
19]. However, in many practical situations, the uncertainties have to be modeled as a larger number of scenarios, which means that the computational time would be unaffordable by iteratively running traditional methods.
Motivated by the requirement of efficiently solving the β-robust p-median problem with a large number of possible scenarios, a novel multi-scenario cooperative evolutionary (MSCE) metaheuristic was proposed. The definition of the neighborhood scenario was introduced, which was used to describe scenarios similar to each other. The major contribution of the MSCE algorithm is to utilize the knowledge obtained in the solutions of neighborhood scenarios. As all of the scenarios are used to model the uncertainties in a p-median problem, there are some similarities between these scenarios. The basic observation is that the optimal solution of a scenario is close to another scenario, if these two scenarios have higher similarities. Thus, the knowledge obtained in the solving of a scenario is helpful to the solving of its neighborhood scenarios. There are two advantages of the MSCE algorithm compared to other approximate heuristics. First, all scenarios are solved parallelly in one run of the MSCE algorithm, which improves the efficiency of computational time. Second, in the parallel search process of MSCE, helpful knowledge is exchanged among neighborhood scenarios for a certain number of iterations, which improves the performance in both optimality and efficiency.
The rest of this paper is organized as follows. The model methodology is presented in
Section 2.
Section 3 proposes the MSCE algorithm and
Section 4 presents the computational experiments. Finally, the paper is concluded in
Section 5.
2. Model Methodology
In this section, a 0–1 integer programming formulation for the β-robust p-median problem is presented. The objective of the problem is to locate the number of p facilities required to minimize the overall expected transportation cost for satisfying all of the demand, subject to the constraints that the relative regret under each scenario is no more than β. In order to calculate the relative regret of a scenario, the optimal solution under this scenario should first be obtained through solving a corresponding deterministic p-median problem. The detail formulation is presented in the following subsections.
2.1. Notations
The main parameters and variables used in the model development are presented as follows.
Indices: |
M | Set of customer nodes, indexed by i; |
N | Set of potential facilities, indexed by j; |
S | Set of scenarios, indexed by s; |
Parameters: |
| the cost of supplying one unit to customer node i from facility j; |
| the demand of customer point i under scenario s; |
p | the number of facilities that can be opened; |
| the probability that scenario s may occur; |
β | a constant that represents the maximum relative regret permitted for each scenario; |
| the optimal objective function value of the problem under scenario s; |
Decision variables: |
| is 1, if facility j is opened in the robust optimization problem; otherwise 0; |
| is 1, if facility j is opened in the problem under scenario s; otherwise 0; |
| is 1, if the demand of customer node i is supplied by facility j under scenario s; otherwise 0. |
2.2. Subproblem Formulation for a Given Scenario
Under a given scenario, the problem is a classical
p-median facility location problem, which can be formulated as follows.
Objective function (1) is to minimize the total transportation cost. In Equation (2), the number of opened facilities under scenario s is restricted to p. Constraint (3) ensures that the demand of the customer node i can only be supplied by opened facilities. Constraint (4) ensures that each customer node must be allocated to a facility. Constraint (5) is a standard integrality constraint.
2.3. Problem Formulation for the β-Robust Situation with Multiple Scenarios
Through solving the model (Ps), we can obtain the optimal solution for the p-median facility problem under any given scenario s, and calculate the corresponding optimal function value, . Given a feasible solution for all scenarios, noted as X, the β-robust solution can be defined below.
Definition 1. Solution X is β-robust for all scenarios when the following constraints are satisfied:whererepresents the objective value for solution X under scenario s. Then, the
β-robust
p-median facility location problem with uncertain demand can be formulated as follows.
Objective function (7) is to minimize the expected overall transportation cost. Submitting function (1) into (6), we can obtain Constraint (8), which ensures that the relative regret for any scenario is no more than the required level β. The meaning of Constraints (9)–(12) is similar to the ones in Ps, while all these constraints must be satisfied for all scenarios .
3. Multi-Scenario Cooperative Evolutionary Algorithm
From the above model methodology, we can see that all of the subproblems
Ps (
) should be solved first, then the problem
β-P can be solved based on the optimal solutions of
Ps. Subproblem
Ps is a typical
p-median problem, which is usually solved by approximate metaheuristics [
1]. When the number of scenarios is huge, it becomes quite time-consuming to obtain all of their (near) optimal solutions by repeatedly running a metaheuristic in the current literature. Thus, a multi-scenario cooperative evolutionary (MSCE) algorithm was developed to calculate the solutions for all scenarios in one run. The basic idea is that the optimal solutions to a set of (
Ps) subproblems have higher similarities if the corresponding scenarios are similar to each other. The scenarios with higher similarities are called neighborhood scenarios. In the MSCE algorithm, the search processes for solving all (
Ps) subproblems are parallelly conducted, and the knowledge of the solutions for similar scenarios is extracted and exchanged among these scenarios for a very certain number of iterations. The utilization of knowledge from the solutions of similar scenarios accelerates the convergence of the algorithm and helps find better solutions.
3.1. Solution Encoding
The solution of the
p-median problem is encoded as a set with
p elements, where each element represents the index of a selected facility for opening. For example, in a 5-median problem with 15 potential facility locations, {2,5,6,12,15} is a code that represents a feasible solution, and the potential locations, 2,5,6,12, and 15, are selected to establish facilities. This encoding strategy has been widely used in GA based algorithms [
4,
9], and can ensure that Constraint (2)/(9) is automatically satisfied.
3.2. Greedy Heuristics for Generating Initial Solutions
In this section, two greedy heuristics are presented to help generate near-optimal initial solutions for the subproblems under each scenario and the β-robust problem.
3.2.1. Greedy Adding Heuristic
The main idea of the greedy adding heuristic (GAH) is to find the potential facility locations sequentially with an attempt that can decrease the most value of the objective function (1), and add them into a set of open facilities one by one until the number of open facilities reaches
p. The idea of the GAH has been utilized to solve different facility location problems [
22,
23].
The main procedure of the GAH is as follows.
1 | Input: the parameters in model Ps; |
2 | the stopping criteria; |
3 | Output: the solution of Ps. |
4 | Step 0: Initialize |
5 | let be the set of opened facilities and initialize ; |
6 | let be the set of unopened facilities and initialize . |
7 | Step 1: calculate the objective value |
8 | For Do |
9 | let and calculate function (1) to obtain for ; |
10 | End for |
11 | Step 2: find and add |
12 | find ; |
13 | add into , and ; |
14 | remove from , and ; |
15 | Step 3: Check stopping criteria |
16 | If the number of elements in < p Do |
17 | go to Step 1; |
18 | Else |
19 | stop the algorithm and output the solution; |
20 | End if |
3.2.2. Greedy Deleting Heuristic
The greedy deleting heuristic (GDH) can be viewed as a reverse process of GAH. First, all potential facilities in set
N are assumed to be opened. Then, the facilities are closed one by one until the number of opened facilities is equal to
p. In each process of closing one facility, the one whose closing can make the minimum objective value is found and deleted from the set of opened facilities. This greedy deleting strategy was also utilized by Berman et al. [
24] and Alp et al. [
4] for solving facility location problems. When the GAH is utilized to solve the
β-P problem, only the one satisfying Constraint (8) can be deleted in each deleting operation.
The main procedure of the GDH is as follows.
1 | Input: the parameters in model β-P; |
2 | the stopping criteria; |
3 | Output: the solution of β-P. |
4 | Step 0: Initialize |
5 | let be the set of opened facilities and initialize ; |
6 | let be the set of unopened facilities and initialize . |
7 | Step 1: calculate the feasible objective value |
8 | For Do |
9 | Let ; |
10 | calculate the function (1) for solution and check Constraint (8); |
11 | If solution can satisfy Constraints (8) for all scenarios Do |
12 | calculate function (7) for solution and note as ; |
13 | End if |
14 | End for |
15 | Step 2: find and delete |
16 | find ; |
17 | delete from , and ; |
18 | Step 3: Check stopping criteria |
19 | If the number of elements in > p Do |
20 | go to Step 1; |
21 | Else |
22 | stop the algorithm and output the solution; |
23 | End if |
3.3. Knowledge Exchange with Neighborhood Scenario
The neighborhood relations between the scenarios are defined based on the distances between their demand vectors.
Definition 2. The distance between two scenarios s andis defined as follows:where s andare the index of scenarios. It can be seen that the smaller the value of , scenario s is more similar to scenario . The optimal solution of the subproblem under scenario s should be similar to that of the subproblem under scenario if scenarios s and are similar to each other. Thus, the knowledge obtained from the solution of subproblem s should be helpful in optimizing subproblem . Based on this observation, three neighborhood knowledge utilizing (NKU) operators are employed to utilize the helpful knowledge in the solutions of the neighborhood scenarios.
(1) NKU operator based on random exchange
Let be the set of opened facilities in a local optimal solution found in the search process of solving subproblem corresponding scenario s and be that of scenario . Let be the common facilities appearing in and , that is . As there are similarities in the optimal solutions of scenario s and its neighborhood scenarios , the common part is thought to be useful and is kept in the new solution. The other part of the solution , that is , may also be part of the optimal solution of scenario s. Thus, some of the uncommon parts in and are randomly exchanged to form new solutions. To illustrate, we considered an instance with 15 potential facilities and p = 5. Let be a local optimal solution of scenario s and be a local optimal solution of its neighborhood scenario . Then , , and . Some of the uncommon parts in and are randomly selected (e.g., 10 in and 8 in ) and exchanged to form two new solutions, and . Finally, the better one of the two new solutions is selected.
(2) NKU operator based on the GAH
In the NKU operator based on the GAH, the common parts of the solutions for scenario s and its neighborhood scenario are also kept in the new solution. The number of facilities in is usually less than p. In this situation, the number of facilities are selected by the GAH from all of the unopened facilities () to form a new solution.
(3) NKU operator based on the GDH
In the NKU operator based on the GDH, the opened facilities in and are united to form a set . Then, the facilities in are deleted one by one using GDH until the number of facilities is equal to p.
The above three NKU operators can exchange some useful knowledge between the solutions of a scenario and its neighborhood. In the MSCE algorithm, these NKU operators are randomly employed.
3.4. Local Search
After the knowledge is exchanged between the solutions of a scenario and its neighborhood, a new solution is generated for the scenario. Then, using this new solution as a starting node, some local searches are conducted to improve the solution. Two local search operators are employed in the local search.
(1) Random exchange
In the random exchange local search operator, a random number of opened facilities in the current solution is randomly selected and exchanged with the same number of facilities randomly selected from the unopened facilities.
(2) Neighborhood exchange
The neighbor exchange local search operator is designed based on a definition of neighborhood facility, which is introduced as follows.
Definition 3. The kth neighborhood facility of facility j is the facility that is the kth closest to facility j among all of the other potential facilities.
In this operator, the randomly selected facilities only exchange open decisions with their neighborhood facilities. For example, are the first five closest facilities of 3. If facility 3 in the solution is selected and exchanged with an unopened facility, one of the unopened facilities in set will be selected for exchange with 3.
3.5. Framework of MSCE
For conducting the MSCE algorithm, the neighborhood scenarios/facilities for each scenario/facility are first found in the initialization step. Then, the GAH algorithm is employed to calculate a local optimal solution for all scenarios, which serve as the starting search point. The GDH algorithm is used to calculate a local optimal solution for the β-P problem, as it has more opportunities to generate a feasible solution that satisfies Constraint (8). In the second step, the useful knowledge on solving a scenario subproblem is extracted and utilized to generate a new solution by randomly employing one of the three NKU operators. In the third step, a number of local search operations are conducted from the new solution obtained in Step 2. The local search process is parallelly conducted and the number of iterations is set as 100. In the fourth step, the local optimal solution obtained in Step 3 is checked to see if it can improve the objective of its neighborhood scenarios. The last step is to check if the stopping criteria are satisfied.
The main procedure of the MSCE algorithm is presented as follows:
1 | Input: the parameters in model Ps and β-P; |
2 | the stopping criteria; |
3 | Ns, the number of neighborhood scenarios. |
4 | Output: the optimal solutions for all Ps () and β-P. |
5 | Step 0: Initialize |
6 | calculate the distances between each two scenarios and obtain their neighborhoods; |
7 | calculate the distances between each two facilities and obtain their neighborhoods. |
8 | Step 1: Generate initial solutions by greedy heuristics |
9 | For i = 1: Do |
10 | calculate the solution for subproblem s by GAH heuristic; |
11 | End for |
12 | calculate the solution for β-P by GDH heuristic. |
13 | Step 2: Exchange knowledge between neighborhood scenarios |
14 | For s = 1: Do |
15 | randomly select a NKU operator to generate a new solution; |
16 | End for |
17 | generate a new solution for β-P by the NKU operator based on GDH; |
18 | Step 3: Parallel local search |
19 | Parfor all Ps () and β-P Do |
20 | For i = 1: T Do |
21 | generate a new solution by a randomly selected local search operator; |
22 | update the solution if the newly generated one is better; |
23 | End for |
24 | End parfor |
25 | Step 4: Update solutions for neighborhood scenarios |
26 | For all Ps () and β-P Do |
27 | For i = 1: Ns Do |
28 | update neighbor i’s solution if the solution of scenario s is better; |
29 | End for |
30 | End for |
31 | Step 5: Check stopping criteria |
32 | If one of the stopping criteria is satisfied Do |
33 | stop the algorithm and output the solutions; |
34 | Else |
35 | go to Step 2; |
36 | End if |
In this paper, the number of neighborhood scenarios, Ns, was set as 20. The number of local search iterations, T, was set as 100. Two stopping criteria were applied: (1) the number of iterations conducting Step 2 was over 1000, and (2) if there was no improvement in the solution for 10 iterations conducting Step 2.
5. Discussion and Conclusions
The main motivation behind the use of the MSCE algorithm was to sufficiently utilize the knowledge obtained from the solutions of neighborhood scenarios to improve the overall search efficiency of solving all scenarios. The computational results for the two practical cases (H-95 and C-668) with 100 scenarios confirmed that the proposed algorithm could obtain better solutions in a much shorter time when compared to the LAH algorithm and GA. The algorithm was suitable for solving problems with a large number of scenarios, as more similar scenarios would be found in a large set of scenarios. If there are a small number of scenarios that are quite different to each other in the problem, there would not be many similarities between these optimal solutions, that is, the MSCE algorithm cannot find sufficient useful knowledge. In this situation, the traditional approximate metaheuristics may be more suitable.
In this work, six instances with different sizes generated from two practical cases were tested and compared in the experiment. In future research, more experiments for different cases should be conducted to further test the performance of the proposed algorithm. The proposed algorithm can also be extended to solve many other optimization problems under uncertain scenarios, for example, the path planning problem in [
27], the portfolio optimization problem in [
28], and the vehicle routing problem with uncertain demand [
29]. Aside from the LAH and GA, other metaheuristics can be compared with the MSCE algorithm. Three neighborhood knowledge utilizing operators are designed to draw and use the helpful knowledge in solutions of neighborhood scenarios. It would be a valuable direction to study new learning strategies for utilizing the similarities among scenarios, then, a more efficient algorithm can be developed.