Formulation and solution of a two-stage capacitated facility location problem with multilevel capacities

In this paper, the multi-product facility location problem in a two-stage supply chain is investigated. In this problem, the locations of depots (distribution centres) need to be determined along with their corresponding capacities. Moreover, the product flows from the plants to depots and onto customers must also be optimised. Here, plants have a production limit whereas potential depots have several possible capacity levels to choose from, which are defined as multilevel capacities. Plants must serve customer demands via depots. Two integer linear programming (ILP) models are introduced to solve the problem in order to minimise the fixed costs of opening depots and transportation costs. In the first model, the depot capacity is based on the maximum number of each product that can be stored whereas in the second one, the capacity is determined by the size (volume) of the depot. For large problems, the models are very difficult to solve using an exact method. Therefore, a matheuristic approach based on an aggregation approach and an exact method (ILP) is proposed in order to solve such problems. The methods are assessed using randomly generated data sets and existing data sets taken from the literature. The solutions obtained from the computational study confirm the effectiveness of the proposed matheuristic approach which outperforms the exact method. In addition, a case study arising from the wind energy sector in the UK is presented.


Introduction
The two-stage facility location problem (TSFLP) with two types of facilities can be classified as a type of hierarchical facility location problem. In the first stage, products B Chandra Ade Irawan chandra.irawan@nottingham.edu.cn 1 produced/supplied by plants are transferred to capacitated depots. The location and the number of plants and depots can be treated as fixed or as decision variables. Their capacity may be finite (capacitated) or unlimited (uncapacitated). In the second stage, the products are delivered to customers. The problem to be addressed includes finding an optimal distribution structure in order to minimise both the fixed (opening) cost of the plants and depots and the transportation costs associated with both stages.
The two-stage location problem has been investigated in the literature. A dual-based optimization procedure for the two-echelon uncapacitated facility location problem was proposed by Gao and Robinson (1992). The two-stage facility location with a single sourcing constraint on depot-plant assignment and customer-depot assignment was investigated by Tragantalerngsak et al. (1997) where six Lagrangean relaxation heuristics are introduced. Marín and Pelegrín (1999) applied Lagrangian relaxation to the resolution of two-stage location problems. Klose (1999Klose ( , 2000 studied the two-stage facility location with a single product, depot location, plant-depot multiple source flow and single source customer-depot assignment. An effective linear programming approach and a Lagrangean relax-and-cut algorithm are proposed to achieve lower and upper bounds for the problem. Tragantalerngsak et al. (2000) also proposed a Lagrangean-based branch-and-bound method to solve the problem. Hinojosa et al. (2000) studied a heuristic algorithm based on Lagrangean relaxation to solve a multi-period two-echelon multicommodity capacitated plant location problem. Keskin and Üster (2007a, b) proposed a scatter search for a multi-type transhipment point location problem with multi-commodity flow and studied meta-heuristic approaches with memory and evolution for a multi-product production/distribution system design problem respectively. Li et al. (2011) proposed a Lagrangean-based heuristic for a two-stage facility location problem with handling costs with multiple products and three layers of nodes: plants with limited production capacities, capacitated depots to be located and customers with known demands per product. The aim of their model is to minimize a total cost comprising depot opening, transportation and handling costs. Li et al. (2014) investigated a multi-product facility location problem in a two-stage supply chain in which plants have a production limit, potential depots have limited storage capacity and customer demands must be satisfied by plants via depots. A hybrid method is developed where the initial lower and upper bounds are obtained by a Lagrangean based heuristic and a weighted Dantzig-Wolfe decomposition and path-relinking combined method are applied to improve obtained bounds. Several variants of the two-stage location problem were also studied by Li et al. (2012), Rodríguez et al. (2014), Camacho-Vallejo et al. (2015), and Mišković and Stanimirović (2016).
The papers cited above deal with the two-stage facility location problem with fixed capacity for each potential depot. However, in practical situations, the capacity of the depot is also considered as a decision variable which needs to be determined (Correia and Captivo 2003). This means that the problem is not only to find the optimal location of the depot but also its capacity. In this study, we not only deal with one product but also with multiple products meaning that a depot may have a different capacity for each product. To the best of our knowledge, this type of problem has not yet been addressed in the literature. Therefore, this paper proposes new mathematical models and a solution method to deal with the two-stage capacitated facility location problem in the presence of multilevel capacities.
The main contributions of this paper are as follows: • Propose for the first time mathematical models for the two-stage capacitated facility location problem in the presence of multi-product and multilevel capacities, • Propose an effective matheuristic approach based on an aggregation method to solve the problem, • Provide a new dataset for the new problem and produce good quality solutions for benchmarking purposes.
The remainder of this paper is organised as follows. Mathematical models for the two-stage capacitated facility location problem considering the presence of multilevel capacities are presented in Sect. 2. Section 3 discusses the proposed matheuristic approach to solve the problem. The computational results are presented in Sect. 4 followed by conclusions in the final section.

Problem formulation
In this section, two mathematical models of the two-stage capacitated facility location problem considering the presence of multilevel capacities are presented. Here, a set of potential depots to choose from is given where the depots to be opened (opened depots) can be determined by solving the models. In the first model, which we refer to as Model A, each potential depot has an associated set of possible capacities for storing each product with different fixed costs. The capacity is related to the maximum amount of units that can be stored for each product in the depot. In this paper, we refer to this capacity as the 'product capacity'. For example, suppose that there are 2 products (Product P1 and P2) where the possible capacities of a depot for these products are 10, 15 and 30 for Product P1 and 80, 120 and 160 for product P2. Here, the first decision is to determine whether we will open this depot or not. The second is to decide whether both products will be stored in this depot. The depot may store both products P1 and P2 or just one of them. Finally, the optimal capacity for each product for this depot needs to be found.
In the second proposed model, termed Model B, the capacity of a depot is based on the size (volume) of the depot that is required to be built. In this model, a set of possible capacities (volume) for each potential depot is given. In this study, we refer to this capacity as the 'volume capacity'. The first decision generated by this model is to determine whether the depot should be opened or not whilst the second is to decide how big a depot needs to be built. Here, we assume that the volume needed to store one unit product is known. The total volume needed to store all products must not exceed the size of the opened depot.

Model A
In model A, there are two types of fixed (opening) cost where the first is the setup (fixed) cost for opening a depot. The second fixed cost is related to the capacity of each product used in the depot. The fixed cost is dependent on the product capacity and the location of the depot. Therefore, the fixed cost of a potential depot may be different from that of others. An opened depot is also not necessarily built to store all products. In other words, the opened depot may keep only selected products. In this model the first total fixed cost can be determined using decision variable Q j and the second one byŶ j pd . The following notations are used to describe the sets, parameters, and decision variables of Model A. Sets I : set of plants with i as its index and l = |I | J : set of potential depots with j as its index and m = |J | K : set of customers with k as its index and n = |K | P: set of products with p as its index and o = |P| D j p : set of product capacities at potential depot j for storing product p with d as its index and δ j p = D j p .
Parameters s i p : the capacity of plant i(i ∈ I ) to produce product p( p ∈ P) w kp : the demand of customer k(k ∈ K ) for product p( p ∈ P) f j : the fixed cost for opening depot j ( j ∈ J ) f j pd : the fixed cost to store product p( p ∈ P) using product capacity d(d ∈ D j p ) in depot j ( j ∈ J ) b j pd : the number of product p( p ∈ P) that can be stored in depot j ( j ∈ J ) when using product capacity d(d ∈ D j p ) c i j p : unit transportation cost of product p( p ∈ P) from plant i(i ∈ I ) to depot j ( j ∈ J ) c jkp : unit transportation cost of product p( p ∈ P) from depot j ( j ∈ J ) to customer k(k ∈ K ) Decision variables The problem can be modelled as an integer linear problem (ILP) as follows.
X i j p ≥ 0, integer ∀i ∈ I, j ∈ J, p ∈ P (7) In the objective function (1), the first term represents the fixed cost of opening the depots, the second term is the fixed cost of the depots to store the products, the third term is the total transportation cost from the plants to the depots and the fourth term is the total transportation cost from the depots to the customers. Constraints (2) ensure that the total number of products transferred from a supplier does not exceed its capacity. Constraints (3) guarantee that the capacity constraints at the depots are satisfied. Constraints (4) indicate that each opened depot only uses at most one capacity level for each product. Constraints (5) ensure that the demand of each customer for each product is met. Constraints (6) state flow conservation constraints for the depots. Constraints (7) and (8) impose non-negativity and integer conditions on the number of products delivered. Constraints (9) and (10) refer to the binary nature of the variables Y and Q (the decisions whether a depot is opened or not and which capacity is used by the opened depot).

Model B
In the second model, model B, the capacity considered is based on the required size (volume) of the depot. In contrast to Model A, the fixed cost in Model B only consists of one term as the fixed cost includes those of both its opening and the storing of products, based on its capacity. In this model the total fixed cost can hence be calculated based on variable decision Y jd . Several possible volume capacities for each depot are considered in this model where the volume capacity of a potential depot has a fixed cost that is also dependent on its size and location. The dimension/volume of each product is required in this model in order to determine the capacity constraints of the depots. The notations used for sets and parameters in this model are similar to the ones provided in the previous model (model A) with some revisions described as follows. Set D j p is replaced byD j whereas parametersf j pd andb j pd are substituted by f jd and b jd respectively. Parameterf j is not required in this model but parameter α p is added. Setŝ D j : set of feasible volume capacities at potential depot j with d as its index andδ j = D j .
Parameters f jd : the fixed cost for opening depot j ( j ∈ J ) using volume capacity d(d ∈D j ) b jd : the volume (size) of depot j ( j ∈ J ) using volume capacity d(d ∈D j ) α p : the volume required to store a unit of product p( p ∈ P) Decision variables X i j p andX jkp as defined in the previous model.
The problem can be modelled as an integer linear problem as follows. Min Subject to j∈J X i j p ≤ s i p , ∀i ∈ I, p ∈ P i∈I p∈P X i j p ≥ 0, integer ∀i ∈ I, j ∈ J, p ∈ P (17) The objective function (11) aims to minimise the sum of the fixed costs of opening depots and transportation costs. The first term of this objective function represents the fixed cost of opening the depots based on the capacity of the depots, the second term is the total transportation cost from the plants to the depots, and the third term is the total transportation cost from the depots to the customers. Constraints (12) enforce the capacity constraints of the plants. Constraints (13) ensure that the size (volume) of depots is enough to store the products. Constraints (14) make sure that each opened depot only uses one volume capacity. Constraints (15) guarantee that the demand of each customer is satisfied. Constraints (16) state flow conservation constraints for the depots.

The solution method
The classical two-stage capacitated facility location problem (TSCFLP) is an NP-hard optimization problem as it represents a generalization of the simple plant location problem, which is proved to be NP-hard by Krarup and Pruzan (1983). The proposed model with multilevel capacities is even harder to solve than the classical TSCFLP using an exact method (via an optimizer software such as CPLEX, Lindo, and Xpress) especially when the size of the problem is relatively large. To overcome this weakness a matheuristic approach is developed by integrating an aggregation technique and an exact method. We refer to this method as a Matheuristic Approach incorporating an Aggregation Technique (MAAT). When the location problems involve a large number of demand points, it may be sometimes impossible and time consuming to solve to optimality (Francis et al. 2009). It is quite common to aggregate demand points/depots when solving large scale location problems. The main idea behind the aggregation is to simplify the problem by reducing the number of demand points/depots to be small enough that an optimiser can be used to solve the reduced problem within a reasonable amount of computing time. However, the approximation involved may lead to a level of sub-optimality when the aggregated solution is put into practice in the actual real-world situation. The aggregation technique has successfully addressed large facility location problems such as for large p-median (Irawan et al. 2014;Irawan and Salhi 2015a) and p-centre problems (Irawan et al. 2016). A review on aggregation techniques for large facility location problems is provided by Irawan and Salhi (2015b). The proposed matheuristic approach (MAAT) is developed to solve both Models A and B presented in Sect. 2. Matheuristics have been successfully used to solve tackle facility location problems (Stefanello et al. 2015;Irawan et al. 2017a, b). The proposed method consists of three stages where the main steps of this approach are depicted in Fig. 1. The first stage is an iterative process that incorporates the aggregation of potential depot sites and the implementation of the proposed local search. Firstly, m potential depot sites are aggregated into μ potential sites, with μ << m. The value of mis determined based on the maximum number (upper bound) of the facilities that need to be opened (ρ). The value of ρ can be approximated by following expressions: Here, μ = β · ρ where β is a parameter. When choosing the aggregated potential depot sites, the aggregation includes the depot sites obtained from the previous iteration (the best solution, S * ). The remaining (μ − |S * |) potential depots are randomly chosen from the mpotential depot sites. The main idea behind this is to make sure that the reduced problem has a feasible solution. The resulting aggregated problem with l plants, μdepots and n customers is then solved by CPLEX within τ seconds. A duality gap (%Gap) is also set as a termination criterion where CPLEX will stop when the %Gap reaches ε%. Let Z be the terminating objective function value and S be the corresponding vector of the obtained facility configuration. The description of the proposed local search is presented in the following subsection. The obtained depots location configuration, if it is better than the previous one, is then fed to the next iteration as part of the set of the aggregated potential depot sites. The process is repeated T times and the best solution (S * ) from this step will be fed to the next step. The values of β and T influence the quality of the solution obtained. The chance of getting a better solution is higher when the values of β and T are set higher as this will increase diversification. However, the computational time also increases for higher values of β and T . In the second stage (Stage 2), the proposed local search is applied to solve the original problem (without aggregation) starting from the best depot configuration obtained from the previous stage. The obtained solution (S * ) from the local search on the original problem is then fed into the final stage where the mathematical formulation of Model A or Model B is solved by an exact method (CPLEX). In the final stage, when solving Model A and Model B, the number of potential depots is reduced from m to |S * |. In other words, the set of potential depot sites (J ) is replaced by S * (the incumbent solution). CPLEX will find the best location to open the depots (if necessary), determine the best capacity for each opened depot, the products flows (integer) from plants to depots and from depots to customers, and calculate the objective function value (Z * ).

The proposed local search
The proposed local search is a hybridisation of the fast interchange heuristic proposed by Whitaker (1983) and an exact method. We enhance the heuristic by incorporating the exact method to solve the multi-product capacitated transhipment problem (MPTP). The exact method is integrated within the local search to optimally solve the transhipment problem whenever the locations of opened depots along with their capacity are known/fixed. Moreover, we also enhance this heuristic by replacing a depot in the current solution with the potential depot (not in current solution) that is not too far from the removed depot. By restricting the search, the local search runs relatively fast at the expense of a small loss in quality.
For Model A, in the case where the location of the opened depots along with their capacity for each product are known, the problem can be treated as the multi-product capacitated transhipment problem, which we refer to as the MPTP-A. The MPTP-A is also relatively easy to solve when we relax the amount of products transported from one node to others to a real value instead of an integer one. The MPTP will hence be a linear programming formulation. Let S ⊂ J be the set of opened depots and a j p be the product capacity used by depot j to store product p. The mathematical model for the MPTP-A is as follows: Decision variables The MPTP-A can be modelled as a linear problem as follows.
For Model B, similarly to Model A, when the location of the opened depots along with their capacities are fixed, the problem can also be treated as the multi-product capacitated transhipment problem which we refer to as the MPTP-B. S is also denoted as the set of opened depots withâ j be the volume capacity used by the depot j. The mathematical model for the MPTP-B is relatively similar to the one for the MPTP-A with minor revisions in objective function (22) and constraints (24). In the MPTP-B, the objective function (22) is replaced by objective function (29) as follows: and constraints (24) are replaced by constraints (30) as follows: i∈I p∈P The main steps of the proposed local search are given in Fig. 2 which is based on the fast interchange heuristic using a first improvement strategy (the exchange process is conducted once there is an improvement). The main objective of the algorithm is to seek a potential depot to be swapped with a one in the current solution where the swap process will be performed if there is an improvement. In this local search, when solving the transhipment problem, all potential depots (J ) are imposed to use the largest capacity for storing each product (for Model A) or the volume (size) of the depot (for Model B).
In this local search, we firstly set the maximum computational time (τ ) to execute the local search. The transhipment problem (MPTP-A and MPTP-B) using the initial solution (S) obtained from the previous step is then optimally solved to evaluate the quality of the solution (Z ). In Steps 3(ii)-3(iii) of Fig. 2, the algorithm aims to allocate the potential depots (not in the current solution) to their nearest opened depot (incumbent solution). The longest/maximum distance (d s ) between the opened depots and their associated potential depots is then determined. The main idea behind this is to restrict the search (the swapping process) by imposing the condition that the substituted depot location must lie within a certain covering radius (γ ·d s ) from the opened depot that will be removed. This will make the local search runs more efficient (in terms of computing time) as the swap process will be skipped when the substituted depot is relatively far from the opened depot that will be removed. In Step 3(iv)a, potential depot j ( j ∈ J, j / ∈ S) is inserted into the solution whereas opened depot sin the current solution is removed. Then, the transhipment problem (MPTP-A and MPTP-B) using S is solved to optimality. In Step 3(iv)d, the swap will be conducted if there is an improvement. The local search will terminate if there is no improvement after all possible swaps based on the incumbent solution have been completed or if the computing time reachesτ seconds.

Computational study
A set of computational experiments have been carried out to evaluate the performance of the proposed solution method. The proposed method was implemented/coded in C++, .Net 2012 where the IBM ILOG CPLEX version 12.6 Concert Library is used to solve the problems with exact method. The tests were run on a PC with an Intel Core i5 CPU @ 3.20 GHz processor and 8.00 GB of RAM. In the computational experiments, two types of dataset are used. The first dataset (Dataset 1) is randomly generated to evaluate our solution method's ability to solve the two-stage capacitated facility location problem considering the presence of multilevel capacities for both models A and B. The second dataset is constructed based on the datasets from the literature and the wind energy industry in the UK. For the case study on the wind industry, the proposed model will be implemented to find the optimal locations for depots required for storing spare parts to support the operation and maintenance of offshore/onshore wind farms.
To evaluate the performance of the proposed matheuristic approach (MAAT), we compare the solutions obtained by the proposed method with those of the exact method (using IBM ILOG CPLEX version 12.63). As the problem is very hard to solve to optimality, the computational time (CPU) for solving the problem using the exact method (CPLEX) is limited to 3 h so the lower bound (LB) and upper bound (UB) can be attained. The performance of the proposed matheuristic method will be measured by %Gap between the Z value attained by the matheuristic approach and the lower bound (LB) obtained from the exact method. Moreover, the %Gap is also set as a termination criterion where CPLEX will stop when the %Gap reached 0.01%. %Gap is calculated as follows: where Z m refers to the feasible solution cost obtained by either the exact method (UB) or the proposed matheuristic approach. In the matheuristic approach, the parameters are set to the following values: T = 10, β = 1.5, τ = 150 s, ε = 0.5%, τ = 108 s,τ = 0.25 · n · |P| seconds and γ = 2.5. For solving Model B, we set the value of β to 2 for |P| = 10 and to 3 for |P| = 5. Those parameters were selected based on a small preliminary study. This selection yields an acceptable performance with respect to the quality of the solution and the computational effort.

Experiments on the randomly generated data (Dataset 1)
In order to conduct extensive computational experiments, we generate a new dataset which we refer to as Dataset 1. This dataset consists of two instances, namely Instance 1A and Instance 1B. Instance 1A is used to evaluate the performance of our method when solving Model A whereas Instance 1B is for Model B. For Instance 1A, there are 20 problems to solve whereas Instance 1B consists of 15 problems. We set the number of products |P| to 5 or 10. The number of plants (l) is varied between 5 and 25 with an increment of 5 whereas the value of mfrom 50 to 500 with an increment of 50. The number of customers (n) is set to 2m with the demand of each customer randomly generated between 1 and 5 for each product. The location of plants, warehouses and customers are generated randomly using a uniform distribution where 2 (n × n) and the coordinates values are integer. The capacity of a plant for each product (s i p ) is generated based on the customer demand. It is assumed that there are three possible capacities for each product for each potential depot (|D j p | = |D j | = 3; j = 1, . . . , m; p ∈ P). We also generate the capacities (capacity/size) of each potential depot for each product (b pd and b jd ) and its fixed (opening) cost along with the transportation cost per km per product. Here, we construct the dataset in such way that in a good solution, the total opening cost is close to the total transportation cost.

Computational results on Instance 1A
The computational results on Instance 1A are presented in Tables 1 and 2, which show the computational results using the exact method (CPLEX) and the proposed matheuristic approach respectively. According to the tables, the complexity of the problem increases when the size of the problem increases as shown by the %Gap value. It is worth noting that when the number of products (|P|) is higher, the problem is more difficult to solve. According to the results shown in Table 1, the problems with |P| = 5 and n ≤ 200 were relatively easily solved by the exact method. In these problems, the %Gap between UB and LB obtained is the requested %Gap termination criterion for CPLEX to solve the problem (i.e. a %Gap of 0.01%). In other words, CPLEX terminated before time based termination criterion of 3 h. Using the exact method, the %Gap value produced is relatively very high when n ≥ 700. On average, the exact method yielded a %Gap of 20.64% which is considered as a large value. The exact method also produced the average portion of fixed (opening) costs of 61.98%. The proposed matheuristic method (MAAT) made a significant improvement in producing solutions on Instance 1A as it provides a better %Gap than the exact method. MAAT produced %Gap of 5.71%, an improvement of almost 15% compared to that of the exact method. The use of MAAT also reduced the average portion of the fixed costs to 41.86, 20.12% lower than that obtained by the exact method. Moreover, MAAT required less than a quarter of the computational time required by CPLEX. Figure 3 shows the location of the opened facilities for problem P1-I3 (n = 200 and |P| = 5) where from 100 potential depots, only 9 need to be opened in order to serve 200 customers. These opened depots will receive 5 types of products from 5 plants and will transfer them to the 200 customers. Here, the demand of a customer for each product is randomly generated between 1 and 5 following a uniform distribution. Table 3 presents the capacity of each plant for each product (p1 -p5) used in problem P1-I3 (n = 200 and |P| = 5) where the location of the plants (x and y co-ordinates) is also given. This problem consists of 100 potential depots, each of which has 3 possible capacities to choose from (30, 50, and 70) for each product. In the solution, only 9 depots are selected to be opened in order to minimise the total cost. Their locations and capacity for each product are given in Table 4.  Fig. 3 The location of opened depots for problem P1-I3 (n = 200 and |P| = 5)

Computational results on Instance 1B
Tables 5 and 6 reveal the computational results on Instance 1B using the exact method and the proposed method (MAAT) respectively. In this instance, CPLEX was not able to solve the problem with n ≥ 600 and |P| = 10 due to memory issues. Therefore, this instance only consists of 15 problems instead of 20. Using the exact method, without the problems with n ≥ 600 and |P| = 10, the %Gap value obtained is relatively low as on average, the exact method provided %Gap of 2.82%. The average proportion of fixed costs from the total cost is 38.29% with the total transportation cost contributing the remainder. Similarly to the previous experiments, the proposed matheuristic method (MAAT) performed very well in solving the problems in this instance. The MAAT produced a %Gap of 2.47%, which is better than that obtained by the exact method. Compared to the exact method, the use of the MAAT decreased the average proportion of the fixed cost by 1.6-36.69%. Similarly to previous experiments, the MAAT also required less than a quarter of the computational time required by the exact method. In general, based on the computational experiments on Dataset 1, the proposed matheuristic technique (MAAT) runs much faster than the exact method while yielding smaller %Gaps.

Experiments on datasets from the literature and wind energy industry
The performance of our proposed approach is also assessed on an existing dataset taken from literature and those used in the wind energy industry.

Experiments on the existing dataset
We test our proposed method on data sets from Eskandarpour et al. (2017) originally used for solving a supply chain network design problem. This existing dataset provides the location of plants, depots and customers where the customer demand for 5 products (|P| = 5) is also given. Here, we use Model A to solve this existing dataset as there is no information related to the dimension of each product. The missing information required to solve Model A is generated based on the total demand for each product. We estimate the capacity of plants and depots along with the fixed cost of opening depot based on its capacity in such way that in a good solution, the total transportation cost is close to the total fixed cost. The existing dataset consists of 15 problems where the number of customer (n) is varied between 60 and 300. Therefore, it can be argued that this existing dataset is relatively small and easier to solve by the exact method than the datasets presented in Sect. 4.1. Tables 7 and 8 show the computational results on the existing dataset using the exact method and the proposed method (MAAT) respectively. According to Table 7, CPLEX terminated before the time based termination criterion of 3 h for 11 of 15 problems where CPLEX stopped because the %Gap between UB and LB has reached the termination level of 0.01%. For the other four problems, the %Gap obtained by CPLEX for these problems is very low after the time based termination criterion of 3 h, indicating near-optimality. On average, the exact method produced a relatively small %Gap of 0.08% with the total transportation cost contributing approximately 40% of the total cost. The proposed matheuristic method (MAAT) also performs well in this instance. The MAAT yielded a %Gap of 0.1% within a relatively short computational time.

Experiments on dataset from the wind energy industry
Renewable energy sources have attracted a lot of attention in recent years due to several factors including a surge in the world energy demand, limitation of fossil fuel reserves, fossil fuel price instability and global climate change (Abdmouleh et al. 2015). The UK Government has set a national target for 15% of its total energy consumption to come from renewable sources by 2020, of which it is expected that wind energy will make the largest single contribution to this target (Jones and Wall 2016). A wind farm can be located either onshore or offshore. The development of the offshore wind industry has significantly increased over the past 20 years. One of the reasons for this growth is that a wind turbine at sea generally produces more electricity than that of its onshore equivalent as the average wind speed at sea is higher (Irawan et al. 2017a, b). The operations and maintenance (O&M) cost is one of the largest components of the cost of a wind farm. One way to reduce the costs is to make the maintenance activities more efficient by optimising the logistic system in order to reduce turbine downtime. The logistic system should hence be designed to ensure that the spare parts are available and easy to be access when they are needed. In the wind energy sector, spare parts are complex and expensive, characterized by high procurement costs and low inventory levels (Tracht et al. 2013). However, it is critical to manage and maintain an adequate level of spare parts as inadequate stocks when a part fails may stop electricity generation and lead to substantial losses.
Spare parts supplied by plants are delivered to capacitated depots, and then distributed to Operation & Maintenance bases (O&M bases). Depots are usually located near to or at the O&M base locations. However, as the inventory levels of the spare parts are relatively low, depots may not be opened at all O&M bases. This means that a depot may serve more than one O&M base. Moreover, a depot may not store the same parts as other depots. For an example, depot A may store only blades and bearing generators whereas depot B may manage transformers and yaw motors.
Optimization of the location and capacity of maintenance accommodations for offshore wind farms has been investigated by a few researchers. De Regt (2012) studied the optimal location of offshore maintenance accommodations by solving a 'Weber' problem to minimise the weighted sum of distances to given points. Besnard et al. (2013) introduced an optimisation model to find the optimal location of maintenance accommodations, number of technicians, choice of transfer vessels and the possibility of using a helicopter to service offshore wind farms.
This section presents computational results of the facility location model for finding the optimal number of depots that need to be opened along with their optimal locations in order to support operation and maintenance of offshore/onshore wind farms in the UK. The capacity to store spare parts for each opened depot is also optimised. Here, Model A is most suitable to be implemented for this case study as the dimension of spare parts and the volume of potential depots are difficult to estimate. In this case, warehouses are treated as depots whereas O&M bases act as customers. The model aims to minimise the total cost which comprises shipment costs (with downtime cost) and capital costs incurred by opening depots. A set of possible product capacities is given where each capacity has a different annual fixed cost which may consist of opening depot and inventory costs. The transfer cost of each component from plants to depots comprises shipping and product costs whereas the one from depots to O&M bases considers downtime and shipping costs. The model will also find the optimal number of depots that need to be opened.
It is common in the wind energy industry that a depot is built to store spare parts of one type of wind turbine. In this case study, the type of wind turbines that we study is Vestas V80/90 as the data for this type of turbine is available in the literature. Therefore, we take into account all wind farm sites (offshore and onshore) in the UK that use the Vestas V80/90 wind turbine. Table 9 shows the detailed wind farm data of sites in the UK that use the Vestas V80/90 (www.renewableuk.com) including the West Gabbard offshore site which is currently still under construction. The table presents the location of wind farms along with number of turbines and installed capacity. Moreover, the table also reveals the location of the associated O&M base for each wind farm. These O&M base locations are also treated as potential depot locations. Table 10 shows the detailed information on parts considered in the case study where the cost and failure rate of each part are given. The part information is based on Lindqvist and Lundin (2010). We also assume that all spare parts are supplied by the manufacturer Vestas located in Randers, Denmark whose coordinates are (Lat 56.433127, Lon 10.047057). In other words, the number of plants is set to one (l = 1).
A wind farm, consisting of several wind turbines with the same specification, is treated as a customer. A wind turbine is formed of several parts which are considered as products. The average frequency of failure of the part per year is defined as failure rate. Therefore, the demand of each part per year for a windfarm is then based on the number of turbines in the windfarm and its failure rate. In other words, the demand of each spare part per year for each wind farm is calculated as the product of the failure rate value and the number of turbines installed in a wind farm. This is acceptable as the number of spare parts needed is determined by how often the part breaks down. Table  10 also presents the possible depot capacities for each spare part, generated based on the demand of each product. The model will select the best capacity for each opened depot.
In the experiments, we also assume that the holding cost per year of each spare part is 20% of its cost. This information is used to calculate the fixed cost. The transportation cost for each component is based on the distance and we set the maximum transportation cost to 20% of the component cost. To implement Model A on this wind energy case study, minor revisions of the mathematical model are needed. First, variable X jkp is treated as a real value instead of an integer as the demand of products for each O&M base is calculated based on failure rates. Second, the equalities on Constraints (6) are replaced by inequalities (≥). This dataset can be solved optimally using the exact method (CPLEX) within a relatively short computational time as the problem is relatively small (l = 1, m = 22, n = 22, and |P| = 36). Therefore, the matheuristic approach (MAAT) is not required to solve this instance.
The optimal solution for this problem reveals that only 4 depots are required to open in order to store the spare parts. Three depots are located on the coast, namely Great Yarmouth, Workington and Ramsgate whereas another depot is located at the Braes of Doune inland windfarm site. This solution is acceptable as there are more offshore than onshore wind turbines. The other main advantage of locating a depot at port is its accessibility from the supplier and customer (O&M base). Moreover, the inland depot is located at the onshore windfarm site that has the largest number of wind turbines. It can also be noted that the locations of the depots are scattered across the UK. The objective function value (the total cost) obtained is 11,451,764.05 where the fixed (opening) cost contributes approximately 35% of the total cost. Table 11 shows the depot configuration located in Port Great Yarmouth, which stores all types of spare parts in the optimal solution.   of opening depots and transportation costs. The first model considers the capacity based on the maximum number of products that can be stored whereas in the second one, the capacity is based on the size (volume) of the depot. As large problems are very hard to solve using an exact method, a matheuristic approach, MAAT (Matheuristic Approach incorporating an Aggregation Technique), is introduced to overcome this weakness. The proposed method is evaluated using a randomly generated dataset and datasets taken from literature and the wind energy sector in the UK. According to the computational experiments, the proposed methods ran efficiently, producing a small %Gap within a short computational time.
The models developed in Sect. 2 can be implemented not only for the wind power sector but also for other industries that need depots to support their business. The models can be enhanced to become bi-objective as there is an underlying trade-off between minimising the number of opened depots and minimising the total costs. The compromise programming method [see Irawan et al. (2015) for more detailed information] can be applied to address the trade-off that occurs. The models can also be extended by considering uncertain customer demand. In this case, a technique such as the stochastic programming could be implemented in order to model the problem.