A MODEL AND TWO HEURISTIC METHODS FOR THE MULTI-PRODUCT INVENTORY-LOCATION-ROUTING PROBLEM WITH HETEROGENEOUS FLEET

. The Multi-Product Inventory-Location-Routing Problem with heterogeneous ﬂeet considers a supply chain, which comprises multiple producers, potential distribution centers (DCs) with opening capacity levels and geographically scattered retailers each of which has deterministic demand over a discrete planning horizon. The goal is determining a set of DCs with their capacity levels to open, assigning retailers to the opened DCs, ﬁnding product quantities to be ordered by and distributed from opened DCs and determining the ﬂeet and routes to satisfy the demands of retailers with minimum cost. A mixed-integer linear programming model is proposed to describe the problem, which is strengthened by two valid inequalities. Since the commercial solver can only solve the very small-sized instances within a reasonable time, two heuristic methods are developed. Results show that the proposed valid inequalities are eﬀective and both methods provide important savings in acceptable run times compared to the commercial solver.


1.
Introduction. According to a study by Blanchard [2], in the USA, more than $1 trillion a year is spent on the optimization of supply chains. To make the relevant optimizations, two main processes, namely production and logistics are generally reviewed. The production process involves mainly the production planning, the handling of materials and the reorganization of the facilities. The logistics process involves the facility location, the procurement of raw materials, the distribution of goods to customers and inventory management issues.
According to a 2005 study by MacroSy research and technology organization, the total cost of the logistics for year in the USA in 2002 was nearly $910 billion, which constituted 8.7% of the country's gross domestic product. Of the $910 billion, $577 billion was for transportation costs, and $298 billion was for inventory holding and the operating costs of distribution centers. The supply chains of the organizations should be designed and optimized in a way that reduces these costs.
Three fundamental problems must be addressed when designing and optimizing supply chains: location-allocation, on the strategic level; inventory management, on the tactical level; and vehicle routing problems, on the operational level. Locationallocation problems are related to choosing a subset of facilities (distribution centers and depots etc.) from candidate facilities and allocating retailers/customers to these facilities in order to minimize the total of facility opening and operating costs. Inventory management problems determine the amount of product(s) for supply and distribution, for each period, in order to minimize the inventory holding and stockout costs. Vehicle routing problems are related to finding the best vehicle fleet and routes in order to minimize fleet and distribution costs.
In many of the previous studies on the aforementioned problems, researchers have studied the problems on their own and in pairs. On the one hand, some researchers claim that the problems should not be addressed simultaneously, as strategic decisions are fixed for long periods of time (more than a year), while tactical and operational decisions are changeable within shorter periods of time.
On the other hand, as noted by Shen and Qi [18]; "These three problems of a supply chain are highly related. An effective distribution scheme depends on the locations of facilities. The allocation of retailers to open distribution centers significantly impacts the optimal order sizes placed at the open distribution centers, as well as each center's routing system", most of the researchers including us claim that the problems should be addressed simultaneously in order to find optimal configurations. Javid and Azad [10] realized this recommendation first and named the considered problem set as the "Inventory-Location-Routing Problem" (ILRP). This problem is one of deciding the location of distribution centers while considering routing and inventory policies over a finite planning horizon.
When we look at the application areas of ILRP, we see medicine/food delivery in humanitarian logistics, ammunition, and other supplies distribution in military logistics, supermarket logistics and many of the logistics cases in which leasing distribution centers are common and distribution to the customers are made in different frequencies. In the long run, when the location costs are in the same order of magnitude as operational costs, as in the pharmaceutical industry, the ILRP comes into play. When direct deliveries might lead to suboptimal solutions, as in the case when vehicle capacity exceeds the maximum demand largely, the ILRP should be considered. We must also use ILRP to perform "what if" analysis during the assessment of supply chain design.
Suppose we think about the case of a catastrophe (earthquake, tsunami and others) where humanitarian logistics must be considered to distribute water, medicine, and other supplies using limited resources. In that case, we face ILRP to locate the facilities where these supplies will be stored/distributed from and determine how to perform distribution under capacity constraints. Failing to provide necessary supplies in time during these situations could cause loss of lives or permanent damages to people's health. Another crucial reason to study ILRP is the increasing importance of greener logistics. To decrease the harm given to the environment and decrease the energy consumed, we should open facilities in convenient places, maintain inventories at these facilities in sufficient levels and minimize the road lengths passed during the distribution process. After making necessary optimizations, we reduce the energy consumed for product storage in cold surroundings and the energy consumed for product distribution. Besides, less carbon dioxide/heat is emitted to the atmosphere.
Regarding the studies on ILRP, Table 1 summarizes the main literature review. The table does not cover the studies regarding fuzzy demand or periodic problems. Columns "# of Pe." and "# of Pr." denote the number of periods and the number of products. Liu and Lee [12] and Liu and Lin [11] considered a single-product, single-period, stochastic demand multi-depot location-routing problem. Liu and Lee [12] proposed a two-phase heuristic method to solve the problem. In the first phase, they find the initial solution through route-first, location-allocation second approach and then in the second phase, improve this solution via the heuristic method. On the other hand, Liu and Lin [11] separated the problem into two subproblems, location-allocation and inventory-routing. They proposed a hybrid heuristic method, which improves both sub problems' initial solutions using tabu search and simulated annealing.
Ambrosino and Scutella [1] proposed models representing different scenarios including plants, central depots, regional depots, transit points, clients and big clients relevant to distribution network design. They tried to solve a single-product, single period, deterministic demand ILRP instance including 13 depots, 95 retailers and 17 vehicles using CPLEX 7.0 within 25 h and reached lower bound with a 23.71% gap.
Ma and Davidrajuh [13] considered a single-product, single period, stochastic demand ILRP and proposed a two-stage iterative solution method. In the first stage, location-allocation decisions are found by solving a mixed integer model. In the second stage, they find inventory decisions using probability theory and routing decisions through a genetic algorithm.
Shen and Qi [18] formulated the single-product, single-period, stochastic demand ILRP as a nonlinear model and proposed a method using Lagrangian relaxation embedded in a branch and bound algorithm. Javid and Azad [10] improved Shen and Qi [18]'s study and presented a model to optimize location, allocation, capacity, inventory and routing decisions simultaneously for the first time. They cast the problem as a mixed integer convex program to solve small-sized instances with branch and bound method. They also proposed a method, which is hybridization of tabu search and simulated annealing to solve medium and large problems.
Granada and Silva [6] considered single-product, multi-period, deterministic demand ILRP and proposed a column generation-based method by decomposing the problem into a master problem (location) and an auxiliary problem (inventoryrouting). They solved small-sized problems with the proposed method, but could not solve big-sized problems because of insufficient RAM.
Sajjadi et al. [16] considered a multi-product, single period, deterministic demand ILRP and proposed a two-stage method. In the first stage, they find location-allocation and routing decisions. In the second stage, they take care of inventory management decisions. They solved instances including two products, ten depots, 200 customers nearly within 10 seconds.
Guerrero et al. [7] presented single-product, multi-period, deterministic demand ILRP and proposed two methods. The first method is sequential, but the second one is hybrid, which embeds an exact approach within a heuristic scheme. They achieved significant savings with hybrid heuristic, but they could only solve relatively small instances compared to the other studies.
Nekooghadirli et al. [14] considered bi-objective, multi-product, multi-period, stochastic demand and bi-objective ILRP. They tried to minimize the total cost and the maximum mean time spent on delivery. They proposed four multi-objective meta-heuristic methods.
Zhang et al. [19] proposed a model for single-product, multi-period, deterministic demand ILRP which makes cross-dock depots assumption. They proposed a hybrid method consisting of initialization, intensification and post-optimization stages. They solved instances including 25 depots, 300 customers and seven periods nearly within 40 minutes.
Guerrero et al. [8] considered the same problem studied by Guerrero et al. [7] and proposed a method: the hybridization of column generation, Lagrangian relaxation and local search. The proposed method founded solutions in shorter computational times than Guerrero et al. [7] study at the expense of solution quality. Ghorbani and Jokar [4] studied multi-product, multi-period ILRP that allows backlogging. They proposed a hybrid heuristic algorithm based on simulated annealing and imperialist competitive algorithm.
Hiassat et al. [9] addressed the multi-period ILRP for perishable products. They solved small-sized instances within acceptable times through CPLEX solver and proposed a genetic algorithm for medium and large-sized instances.
Finally, Saragih et al. [17] considered single-product, probabilistic demand, multi-period ILRP in a three-echelon supply chain. They proposed a heuristic method consisting of two stages. In the first stage, the initial solution is constructed and in the second stage, it is improved iteratively using simulated annealing.
Considering the studies on ILRP up to now, we see that ILRP is dealt with different assumptions in different studies. A comprehensive/generic approach, which takes care of many aspects of the problem altogether is needed. Therefore, in this study we consider multi-product ILRP (MILRP) with heterogeneous fleet, which simultaneously takes care of many aspects of the problem. Main contributions of this study are threefold: first, a mixed-integer linear programming model is proposed to formulate the considered problem. Second, a sequential heuristic method and third a hybrid heuristic method is developed to solve all sized instances of the MILRP with heterogeneous fleet.
The outline of the paper is as follows. In the next section, the mathematical model for MILRP with the heterogeneous fleet is presented first, and then two valid inequalities proposed for ILRP are adapted to the MILRP. In the third section, the sequential heuristic method and in the fourth section the hybrid heuristic method is presented in detail. We give the computational results in the fifth section. The last section includes the conclusions and future research.
2. Problem definition. The problem studied in this paper considers a supply chain design consisting of multiple producers, potential distribution centers (DCs) and a set of retailers for a finite-horizon with a limited number of periods. The goal is to determine a set of DCs with their capacity levels to open, assign retailers to the opened DCs and for all periods find the product quantities to be ordered by and distributed from opened DCs and determine the routes to satisfy the demands of retailers with minimum cost. Representation of the problem is given below: The assumptions of the problem are as follows. Each retailer faces a deterministic demand for each type of product for multiple periods. Each producer produces only a single product. Each DC can be opened in one of the capacity levels (small, medium and large) with different costs. Location and allocation decisions are strategic and fixed for all periods. Inventories are kept both at the DCs and at retailers. There can be initial inventories at the DCs and retailers. There is a holding plus obsolescence cost for each type of product kept at DCs and retailers. Backlogging or stock-out is not allowed. In any period, each vehicle travels at most on one route, and retailers are visited at most once. Vehicles depart from and return to the DCs. Product delivery from producers to DCs is dedicated, whereas from DCs to retailers is non-dedicated. Deliveries are made by heterogeneous vehicle fleet, which contains three types (small, medium and large) of vehicles. The distribution cost includes travelling distance related cost and vehicle fixed cost (using a vehicle at least once cost). Each vehicle type has a different fixed cost and unit cost of travelling per km. After solving the problem, we will get the answers to these questions: • Which DCs should be opened? • dist ij d jtn : demand of retailer j for product n at period t (∀j ∈ J, ∀t ∈ H, ∀n ∈ N ) , OC in : ordering cost of product n for DC i (∀i ∈ I, ∀n ∈ N ) , h ktn : total unit holding with obsolescence cost of node k at period t for product n (∀k ∈ K, ∀t ∈ H, ∀n ∈ N ) , OP C e i : opening cost of DC i with capacity level e (∀i ∈ I, ∀e ∈ E) , M CQ i : periodic operating cost quotient of DC i (∀i ∈ I) , W e i : capacity of DC i opened at capacity level e (∀i ∈ I, ∀e ∈ E) , W j : capacity of retailer j (∀j ∈ J) , Q b : capacity of type b vehicle (∀b ∈ B) , F C b : cost of using type b vehicle at least once (∀b ∈ B) , IS kn : initial inventory of node k for product n (∀k ∈ K, ∀n ∈ N ) , α : weight factor for ordering cost, β : weight factor for distribution cost, θ : weight factor for inventory and obsolescence cost total for DCs, γ : weight factor for inventory and obsolescence cost total for retailers.

Decision variables
   1 if the arc (i, j) is crossed from i to j by vehicle v of type b at period t, (∀i, j ∈ K, ∀t ∈ H, ∀b ∈ B, ∀v ∈ V ) 0 otherwise s itn : quantity of product n ordered by DC i at period t (∀i ∈ I, ∀t ∈ H, ∀n ∈ N ) , w b ijtnv : quantity of product n delivered from DC i to retailer j at period t by type b vehicle (∀i ∈ I, ∀j ∈ J, ∀t ∈ H, ∀n ∈ N, ∀b ∈ B, ∀v ∈ V ) , q ktn : quantity of product n in inventory of node k during period t ∀k ∈ K, ∀t ∈ H 0 , ∀n ∈ N , nv b i : maximum number of type b vehicles used by DC i in the planning horizon (∀i ∈ I, ∀b ∈ B) , u b jtv ∈ {0, 1, . . . , |J| − 1} : subtour elimination variable (∀j ∈ J, ∀t ∈ H, ∀b ∈ B, ∀v ∈ V ) , i∈I j∈J j∈J v∈V Regarding objective function, equation (1) is the sum of the opening plus operating cost, using vehicles at least once cost, ordering cost, distribution cost and holding with obsolescence cost. The set of Eq. (2) and (3) are inventory balance equations for DCs and retailers. The set of Eq. (4) and (5) assigns the initial inventory levels to DCs and retailers. We guarantee the capacity of DCs with constraints (6), the capacity of retailers with constraints (7) and the capacity of vehicles with constraints (8). Constraints (9) ensure the ordering decisions at DCs. Constraints (10) state that retailers can be replenished only by their assigned DCs. Constraints (11) and (12) guarantee the vehicle visits for replenished retailers. Each retailer is assigned to a single opened DC by constraints (13) and (14). Constraints (15) state that each DC can be opened in one of the capacity levels.
Regarding the distribution, constraints (16)-(25) ensure feasible routing. Equations (16)- (18) are traditional vehicle flow conservation constraints. Equations (19) ensure that each vehicle can perform at most one route per period. There cannot be any route between DCs as enforced by equations (20). Constraints (21) find the number of vehicles used by DCs. Constraints (22) state that a retailer can be linked to a DC if it is assigned to that DC. Routes must be related to the opened DCs as enforced by constraints (23) and (24). Constraints (25) are for sub tour elimination. Constraints (26)-(34) state the nature of decision variables.
We guarantee the triangle inequality (dist ij ≤ dist ik + dist kj ) through the parameters used in instance creation, otherwise visits without distribution should be forbidden by adding constraints (35) to the model.

2.2.
Valid inequalities. M1 model can be strengthened by additional inequalities. However, inequalities for the traveling salesman problem, vehicle routing problem and inventory-routing problem are not valid as stated in Guerrero et al. [7]. Nevertheless, we adapted the two valid inequalities proposed for ILRP by Guerrero et al. [7] to the MILRP with heterogeneous fleet as follows.
Theorem 2.1. The inequalities (36) are valid for the M1 model.
Proof of Theorem 2.1. If the DC i distributes products to satisfy the demands at period t -i.e., n∈N j∈J b∈B v∈V w b ijtnv > 0 -then: • At least n∈N j∈J b∈B v∈V w b ijtnv /max b∈B Q b vehicles must depart from DC i at period t to satisfy such a demand -i.e., • At least one vehicle must depart from DC i at period t to satisfy such a demand -i.e., b∈B v∈V , if the capacity of the largest vehicle is loose ( min max b∈B Q b , max e∈E (W e i ) = max e∈E (W e i ) ) . Theorem 2.2. The inequalities (37) are valid for the M1 model.
Proof of Theorem 2.2. The set of constraints (37) state that the minimal number of times a retailer must be visited up to period t equals to the total demand of that retailer that cannot be satisfied by its initial inventory divided by the capacity of the largest vehicle.
3. Sequential heuristic. We see that only very small-sized instances can be solved exactly through IBM ILOG CPLEX solver within a reasonable computation time after doing preliminary experiments. Therefore, we develop the sequential heuristic method (SHM) to solve all sized MILRP instances within a reasonable time. Table 2 gives a general overview of the SHM. In the SHM, we find the location-allocation and inventory management decisions deterministically through the f indSeqSolW ithoutRoutes procedure first. Then, we find the routing decisions stochastically through f indAllRoutes procedure. Hence, only the routing decisions differ between runs.
In the f indSeqSolW ithoutRoutes procedure presented in Table 3, we find location-allocation decisions through f indLocAllocDecisions procedure first and then, we find the inventory management decisions through the f indInvM anDecU singInvM anM odel procedure. Finally, we check the capacity levels of opened DCs and set the minimum required capacity level to each opened DC through the checkAndU pdateOpenedDCsCapLevels procedure. Table 2. Pseudocode of the sequential heuristic method Procedure: solveSeq 1: S 00 ← findSeqSolWithoutRoutes() 2: for repNum = 1 to numOfReps do // for all replications 3: S 0 = S 00 4: S * ← findAllRoutes(S 0 ) 5: Save the S * as the best solution found for the replication numbered repNum Table 3. Pseudocode of the procedure for finding sequential solution without routes Procedure: findSeqSolWithoutRoutes 1: S 00 ← findLocAllocDecisions() 2: S 00 ← findInvManDecUsingInvManModel(S 00 ) 3: S 00 ← checkAndUpdateOpenedDCsCapLevels(S 00 ) 4: return S 00 3.1. Finding the location-allocation decisions. In the f indLocAllocDecisions procedure, we first calculate the estimated total costs for each DC to make location decisions. The estimated total cost for a DC is the sum of estimated opening, estimated ordering, and estimated distribution costs for that DC.
For the estimated opening cost for each DC, we take the average of its opening costs. For the estimated ordering cost for a DC named i we do the following calculation. First, we assume that all retailers are assigned to DC i. Then, for each product, we find the number of required orders to satisfy the total demand considering the initial stock of DC i and multiply this number with the ordering cost of product n for this DC. After summing these values, we find the estimated ordering cost for DC i.
The estimated distribution cost for DC i is calculated as follows. For each retailer, we find the number of required distributions times travelling cost to satisfy the retailer demand with minimum capacity vehicle. Then we sum these costs and find the estimated distribution cost for DC i.
After finding the estimated total costs for DCs, DCs are opened by starting from DC with the minimum estimated total cost. While assigning a retailer to a DC, we assume that the retailer's demands are supplied and delivered at the same demand periods, and we check the DC's capacity accordingly. When an opened DC's highest capacity level is exceeded, next DC is opened and filled with the remaining retailers. While assigning the retailers to an opened DC, retailers are sorted in ascending order according to their distances to the opened DC and then assigned in this order. When all the retailers are assigned, the location-allocation process finishes.

3.2.
Finding the inventory management decisions. We developed the M2 inventory management model to find inventory management decisions exactly. Through the f indInvM anDecU singInvM anM odel procedure, location-allocation decisions are given to the M2 model as input. The M2 model is then created and solved to find the product quantities to be ordered and delivered by opened DCs for all periods. An important assumption made in the M2 model is using only the minimum capacity vehicle for distribution to reduce solution time. This assumption is made with no doubt because, as stated in instance creation parameters, the maximum demand of all retailers can be delivered with the minimum capacity vehicle. M2 model is as follows: , c 0 : unit cost of travelling per km for minimum capacity vehicle, c 0 ij : cost of travelling for i − j arc passed by minimum capacity vehicle (∀ (i, j) ∈ A) , c 0 ij = c 0 . dist ij d jtn : demand of retailer j for product n at period t (∀j ∈ J, ∀t ∈ H, ∀n ∈ N ) , OC in : ordering cost of product n for opened DC i (∀i ∈ I o , ∀n ∈ N ), h ktn : total unit holding with obsolescence cost of node k at period t for product n (∀k ∈ K, ∀t ∈ H, ∀n ∈ N ) , W E i : capacity of DC i opened with capacity level E i (∀i ∈ I o , ∀E i ∈ E), W j : capacity of retailer j (∀j ∈ J), Q 0 : capacity of minimum capacity vehicle, F C 0 : cost of using minimum capacity vehicle at least once, IS kn : initial inventory of node k for product n (∀k ∈ K, ∀n ∈ N ) , α : weight factor for ordering cost, β : weight factor for distribution cost, θ : weight factor for inventory and obsolescence cost total for DCs, γ : weight factor for inventory and obsolescence cost total for retailers,

Decision variables
s itn : quantity of product n ordered by opened DC i at period t (∀i ∈ I o , ∀t ∈ H, ∀n ∈ N ) , w Ij jtn : quantity of product n delivered to retailer j from its assigned DC I j at period t ∀I j ∈ I o , ∀j ∈ J Ij , ∀t ∈ H, ∀n ∈ N , q ktn : quantity of product n in inventory of node k during period t ∀k ∈ K, ∀t ∈ H 0 , ∀n ∈ N , nv 0 i : maximum number of minimum capacity vehicles used by opened DC i in the planning horizon (∀i ∈ I o ), Regarding the objective function; equation (38) is the sum of using minimum capacity vehicles at least once cost, ordering cost, distribution cost and holding with obsolescence cost. The set of Eq. Two valid inequalities adapted to the M1 model are also adapted to the M2 model as follows.
Theorem 3.1. The inequalities (55) are valid for the M2 inventory management model.
Proof of Theorem 3.1. If an opened DC i distributes products at period t to satisfy the demand of retailers assigned to it -i.e., n∈N j∈J I j w Ij jtn > 0 -then: • At least n∈N j∈J I j w Ij jtn / Q 0 number of deliveries must be made by DC i at period t to satisfy such a demand -i.e., j∈J I jxIj jt ≥ n∈N j∈J I j w Ij jtn / Q 0 , if the capacity of the minimum capacity vehicle is tight (min Q 0 , W • At least one delivery must be made by DC i at period t to satisfy such a demand -i.e., j∈J I jxIj jt ≥ n∈N j∈J I j w Ij jtn / W Ij , if the capacity of the minimum capacity vehicle is loose min Q 0 , W Proof of Theorem 3.2. The set of constraints (56) state that the minimal number of deliveries made to a retailer up to period t equals to the total demand that cannot be satisfied with the initial inventory of that retailer divided by the capacity of minimum capacity vehicle.
3.3. Finding vehicle routes. Vehicle routes of the sequential solution are found using an algorithm based on Pasha et al. [15]. The significant differences between our algorithm and the referenced algorithm are explained in the relevant parts. The f indAllRoutes procedure presented in Table 4 finds the vehicle routes for all open DCs in the input solution (sol).
In the procedure first, for each open DC, we find the routes for all periods independently. Hence, we obtain routes requiring different (same with a small probability) vehicle fleets for each period. Next, for each open DC we find the fixed fleet (fleet that will be used for all periods) using four pre-defined fixed fleet finding strategies, and if the found fixed fleet has never been tried before, we run the following process. The aim of using these strategies is to vary the tried fixed fleets, balance the fleets required for each period and decrease the distribution plus using vehicles at least once  Assign the next open DC to openDC i 3: for t = 1 to numOfPeriods do // For all periods 4: If there is a distribution (dist.) to openDC i at period t 5: Find the routes for openDC i at period t using findRoutes procedure 6: Add the fleet of openDC i to the new created triedFleetsForOpenDCI list 7: for strInd = 1 to numOfStrategies do //For all of the fleet finding strategies 8: Find a new fleet (newFleetForOpenDC i ) for openDC i with the strategy in turn

9:
If newFleetForOpenDC i is not included in the triedFleetsForOpenDCI list 10: Add newFleetForOpenDC i to the triedFleetsForOpenDCI list 11: for t = 1 to numOfPeriods do // For all periods 12: If there is a distribution to openDC i at period t 13: If the fleet of openDC i at period t is not a subset of newFleetForOpenDC i 14: Find and assign routes to openDC i for period t with newFleetForOpenDC i

15:
If all of the new routes of openDC i is feasible 16: If dist. plus fixed cost of new routes for openDC i is less than the original cost 17: Assign the new routes to openDC i The f indAllRoutes method presented in Table 4 differs from the method proposed by Pasha et al. [15] in steps 6-20. The most crucial difference is using four strategies consecutively in this study instead of using just one randomly generated strategy as in the Pasha et al. [15]. The first three strategies are Average Plus Large Vehicle Capacity (ALVC), Average Plus Highest Frequency Fraction Vehicle Capacity (AHFFVC) and Vehicle Mix from the Largest Day (VMLD) strategies that give the first three best results as stated in Pasha et al. [15]. ALVC and AHFFVC strategies find the average number of each vehicle type over all periods and add vehicles until the required capacity is reached. The VMLD strategy uses the same fleet as found by the best solution on the period that has the highest delivery amount. The algorithmic details of these strategies are given in the referenced study. We propose a fourth strategy which is the rounded-up form of the ALVC strategy.
We find the routes of an open DC for a specific period t under the constraint of fixed fleet or not with f indRoutes procedure presented in Table 5. The procedure is based on the tabu search meta-heuristic proposed by Glover and Laguna [5]. Neighborhood solutions are found through the shift and swap moves. Aspiration criterion in tabu search defines the situations when the tabu criterion can be overridden. In our implementation, we select aspiration criterion as leading a new best solution. Infeasible solutions, in which the amount of product transported on a route is greater than the assigned vehicle's capacity, are also allowed and penalized with the beta parameter. ALF A parameter updates the beta parameter according to the current solution's feasibility.
With SEH (shrinking and expanding) sub-procedures used in the f indRoutes procedure, routes taken as input are shrunken or expanded for diversification purpose. We use the filling degrees of routes to make shrinking or expanding decision. When a route's filling degree is above the pre-determined filling degree threshold (fillingDegreeThres), we first assign higher capacity vehicles consecutively until the filling degree is dropped under the f illingDegreeT hres. After assigning the highest capacity vehicle to the route, if the filling degree does not fall under the f illingDegreeT hres, we divide the route into new routes. On the other hand, when a route's filling degree is below the f illingDegreeT hres, we first assign lower capacity vehicles consecutively until its filling degree exceeds the f illingDegreeT hres. After assigning the lowest capacity vehicle to the route and its filling degree does not exceed the f illingDegreeT hres, we empty the route by assigning all of its retailers to the other nearest route(s). These operations are made under the fixed fleet constraint in SEHW ithF ixedF leet procedure and without the fixed fleet constraint in SEHW ithoutF ixedF leet sub-procedure. SEH subprocedures used in our study are similar to the SEH procedures used in Pasha et al. [15].
The f indRoutes procedure differs from the method proposed by Pasha et al. [15] with the following major differences. The first difference is our implementation's extra capability that it can also be run under the fixed fleet constraint. The second difference is considering four neighborhood solutions per iteration instead of two to search the neighborhood more effectively. The third difference is the permission of infeasible solutions if the vehicle capacity is not exceeded over a pre-defined ratio. With this capability, we rule over the search space and shorten the computation time. The fourth difference is the continual update of tabu tenure and SEH procedure call frequency with respect to the global iteration counter used in the procedure. With this difference, we try to escape from local minimums, search the neighborhood more effectively and get better solutions.
With f indInitialRoutes procedure, we find the initial routing solution of an open DC for a specific period under the fixed fleet or no constraint. If the procedure is called without the fixed fleet constraint, we create a new route and assign a random Table 5. Pseudocode of the procedure that finds the routes for a specific DC at a given period If(f ixedF leet == null) // Apply the appropriate SEH procedure to S 0 18: Else 20: S 0 ← SEHWithFixedFleet(S 0 , t, f ixedF leet) 21: If(isF easible(S 0 ) && (S 0 .getCost() < S * .getCost())) 22: S * = S 0 23: sEHCallCounter = sEHCallCounter + 1 24: If(isF easible(S 0 )) // Update the beta parameter appropriately 25: beta = beta * (1 -ALFA) 26: Else 27: beta = beta * (1 + ALFA) 28: return S * vehicle type to this route first. Then we fill this route with retailers that must be replenished at that period. When the vehicle's capacity is filled, we create a new route, assign a random vehicle type to this route and continue the same process until all retailers who should be replenished at that period are assigned to a route. If the procedure is called under the fixed fleet constraint, vehicles assigned to the routes are selected from the fixed fleet and the same process is run. When all the vehicles in the fixed fleet are used and if any retailers without a route assignment exist, they are assigned to the routes, which are sorted in ascending order according to filling degree. So, the feasibility of the routes under fixed fleet constraint is not guaranteed.
4. Hybrid heuristic. To find better quality solutions than the sequential heuristic method, decrease the solution time and increase the size of the problems that can be solved in acceptable running times, the hybrid heuristic method (HHM) is developed and used through solveHybrid procedure. The pseudocode of the HHM is presented in Table 6. Empty the triedLocAllocConfigsList 4: S 0 = S 00 // Assign the base solution to initial solution 5: Add the location allocation decision of S 0 to the triedLocAllocConfigsList 6: S 0 ← findAllRoutes(S 0 ) // Find all of the routes for initial solution 7: S * ← intensify(S 0 , triedLocAllocConfigsList) 8: S new = S * 9: diverCounter = 0 // Initialize diversification counter 10: while (diverCounter < totalNumOfDiversifications ) do 11: S new ← diversify(S new , triedLocAllocConfigsList) 12: S new ← intensify(S new , triedLocAllocConfigsList) 13: If (S new .getCost() < S * .getCost()) 14: S * = S new 15: diverCounter = diverCounter + 1 16: Save S * as the best solution found for the replication numbered repN um In the solveHybrid method, we create the initial solution without routes through the f indHybridSolW ithoutRoutes procedure and use this solution in all replications. In a replication first, the base solution is assigned to the initial solution. This solution's routes are found through the f indAllRoutes procedure, which is the same procedure used in the sequential method. This solution is then improved through the intensif y procedure and the best solution found is assigned to the solution S new . Finally, we run the loop, which calls diversif y and intensif y procedures sequentially for totalN umOf Diversif ications times and saves the best solution found for the relevant replication. In the loop, we change the location and allocation decisions of the solution S new through diversif y procedure and then improve the new solution through intensif y procedure. In order to create different solutions (having different location and allocation decisions), the location and allocation configurations tried before are stored in the triedLocAllocConf igsList which is used and updated by intensif y and diversif y procedures. In the diversif y procedure a new solution is created until we make sure that it is not tried before (not stored in the triedLocAllocConf igsList).

4.1.
Creation of the Initial Solution. As stated previously, we create the initial solution without its routing decisions through the f indHybridSolW ithoutRoutes procedure. In this procedure first of all, we find the location and allocation decisions of the initial solution by using the same algorithm used in the sequential heuristic. Then, we find the inventory management decisions of the initial solution by letting the quantities of products delivered to the retailers are equal to the corresponding demands of that period if they cannot be supplied from the retailers' initial stocks. Open DCs in the initial solution order products from producers after the consideration of their initial stock levels. Finally, we find the routing decisions of the initial solution through the f indAllRoutes procedure. To find the routing decisions fast, we call the f indAllRoutes procedure with 0.1 parameter, which is multiplied by globalCounter parameter used in this procedure.

4.2.
Intensification. In the intensification process, the intensif y procedure presented in Table 7 is used to improve the given solution. This procedure is the hybridization of the simulated annealing and the tabu search. We divide the intensification process into phases by using the initial and freezing temperatures (INI TEMP and FRE TEMP) of the simulated annealing.
At the beginning of the procedure, we find the phase number, which the procedure is currently in.
Then, we call the setIntP arams subprocedure with this phase number to assign the appropriate values to the probability parameters used in intensification moves. At each temperature, NUM OF ITERS AT ALL TEMPS repetitions are made. If the current solution is not improved consecutively for NUM OF CONS TEMPS SOL NOT IMP times, a random number is generated. If the value of this number is less than or equal to the ASSIGN THE BEST SOL RATIO, the best solution is assigned to the current solution, otherwise ordering decisions of the open DCs are optimized through updateOrdersW ithOrderingM odel sub-procedure.
There are nine intensification moves. The first four of these moves shift deliveries. At the end of the intensification procedure, we first find the vehicle routes of the best solution again.
Then, we optimize the capacity levels and orders of the open DCs through optimizeDCsCapLevelsAndOrders subprocedure. In the optimizeDCsCapLevelsAndOrders sub-procedure, we compare the best solutions found with two different methods and return the better one.
In the first method, we optimize the orders of open DCs through updateOrdersW ithOrderingM odel sub-procedure and then minimize the capacity levels of all open DCs. In the second method, we bring down the capacity levels of open DCs one level (if there is any) and then optimize their orders one by one through the updateOrdersW ithOrderingM odel sub-procedure. If the updateOrdersW ithOrderingM odel sub-procedure finds a feasible solution for the DC taken as input, DC's new capacity level remains, otherwise its capacity level is changed to the original one. 4.2.1. Shift Delivery Move 1. In this move applied by applyShiftDeliMove1 subprocedure, a random retailer j and then a random delivery period t 1 for this retailer is selected. The delivery made at period t 1 is tried to be shifted to the Table 7. Pseudocode of the intensification process If (S m .getCost() < S * .getCost()) // If the best solution is improved 39: S * = S m 40: temp = temp * COOLING SPEED // Update the current temperature 41: findAllRoutes( S * ) // Find all of the routes for the best solution again 42: S * ← optimizeDCsCapLevelsAndOrders(S * ) 43: return S * randomly selected new period t 2 . We apply this move when the move is not tabu regarding the retailer j, t 1 and t 2 periods. In the move, we find the amounts of products for each product type, which will be shifted to the period t 2 first. While finding these amounts we consider if there is a route to the retailer j at t 2 , if the vehicle type of this route could be changed, if a new route specific to the retailer j on t 2 should be created. After these considerations, we also check the capacity of the DC to which retailer j is assigned (I j ), the capacity of the retailer j and stock out situations. Then, we adjust the amounts of products for each type, which will be shifted to the t 2 period. After these, we also adjust the ordering and delivery plan of the DC I j . Finally, we make the necessary updates to the routes for t 1 and t 2 periods. According to the pre-determined probabilities, we bring down the capacity level of the DC I j and find the routes departing from DC I j again. We see that this move's applicability ratio is around 0.4 and the improvement ratio of the best solution after applying this move is around 0.03.

4.2.2.
Shift Delivery Move 2. The move applied by applyShiftDeliMove2 subprocedure is the two-shift delivery operator proposed by Zhang et al. [19]. The representation of this move is given below. In this move, we select two random retailers, j 1 and j 2 which share the same route for a random period t 2 . We then try to shift the delivery made to the retailer j 1 at period t 2 to the period t 3 . After making space in the shared route, we try to shift the delivery made to the retailer j 2 at period t 1 to the period t 2 . As we can see, this move is the consecutive application of shift delivery move 1. We apply these moves if the moves are not tabu regarding the retailers and periods. According to the pre-determined probabilities, we bring down the DC's capacity level to which retailers j 1 and j 2 are assigned and find the routes departing from this DC again.
Regarding the results, we see that this move's applicability ratio is around 0.2 and the best solution's improvement ratio after applying this move is around 0.02.

Shift Delivery Move 3.
The move applied by applyShif tDeliM ove3 subprocedure is the shift-delivery-deepest operator proposed by Zhang et al. [19]. In this move a random retailer j and then a random delivery period t 1 for this retailer is selected. Then, we try to shift the total delivery made at period t 1 to the other delivery periods of this retailer. Before the shift operations, we order the other delivery periods in decreasing the delivery amount. Delivery shift operations are made as in the shift delivery move 1. We apply this move is if the move is not tabu regarding the retailer j and t 1 period. According to the pre-determined probabilities, we bring down the DC's capacity level to which retailer j is assigned and find the routes departing from this DC again.
Regarding the computational results, we see that this move's applicability ratio is around 0.6 and the best solution's improvement ratio after applying this move is around 0.03.

Shift Delivery Move 4.
The move applied by applyShif tDeliM ove4 subprocedure is the consecutive application of the shift delivery move 3 to the retailers that share a randomly selected route. This move aims to empty and delete a selected route. According to the pre-determined probabilities, we bring down the DC's capacity level to which retailers are assigned and find the routes departing from this DC again. Regarding the computational results, we see that this move's applicability ratio is around 0.34 and the best solution's improvement ratio after applying this move is around 0.03.

4.2.5.
Shifting a Retailer Move. When there are more than one open DCs in the current solution, with the applyShif tRetaM ove sub-procedure, we select a retailer j randomly and try to assign this retailer to an open DC I k , which is different from its assigned DC (I j ). In this move, we cancel the orders and deliveries made by DC I j to supply retailer j's demands first. Then, the deliveries made by DC I j to retailer j are tried to be made by DC I k . We try to make the deliveries in the original amounts and at the original delivery periods during this process. If we could not make the deliveries precisely at the original periods, we try to make the remaining deliveries at other periods. While doing this process, we check stock outs, the capacity of retailer j, and the capacity of DC I k . After making the necessary adjustments, it is possible to have an infeasible solution. If we get a feasible solution, we create a new location-allocation configuration for the new solution and look for it in the triedLocAllocConf igsList. While comparing two configurations, we check if they have the same open DC(s) first. If that is, then we compare their assigned retailers. While comparing the assigned retailers, we try to see that half of the assigned retailers of two compared DCs match. If that is so for all open DC(s), then we say that these two configurations are equal. We add the new configurations to the triedLocAllocConf igsList. According to the pre-determined probabilities, we bring down the capacity level of the DC I j and find the routes departing from DCs I j and I k again. Regarding the computational results, we see that this move's applicability ratio is around 0.4 and the best solution's improvement ratio after applying this move is around 0.001. 4.2.6. Swapping Two Retailers Move. When there are more than one open DCs in the current solution, with applySwapRetaM ove sub-procedure two randomly selected retailers j 1 and j 2 , which are assigned to different open DCs I j1 and I j2 are tried to be swapped. This move is the consecutive application of two shift retailer moves, so we apply the same methodology. In this move, after making the necessary adjustments, it is possible to have an infeasible solution. If we get a feasible solution, we create a new location-allocation configuration for the new solution and look for it in the triedLocAllocConf igsList. If the list does not contain the new configuration, we add it to the list. On the other hand, if we get an infeasible solution, we cancel this move's application. According to the pre-determined probabilities, we bring down the capacity levels of the open DCs I j1 and I j2 and find the routes of these open DCs again.
Regarding the computational results, we see that this move's applicability ratio is around 0.25 and the best solution's improvement ratio after applying this move is around 0.0003.

4.2.7.
Finding All of the Routes Move. In this move, which is applied by f indAllRoutes sub-procedure, we find the routes of the open DCs again. Regarding the computational results; we see that the improvement ratio of the best solution after applying this move is around 0.002.

4.2.8.
Changing the Open DCs' Capacity Levels Move. This move, which is called by changeDCsCapLevelsIf N ecessary sub-procedure, is in fact a combination of moves and could only be applied when the intensification phase has changed. In the move, the capacity levels of open DCs whose filling ratio is less than or equal to the DC FIL DEGREE MIN LEVEL are tried to be brought down and the capacity levels of the opened DCs whose filling ratio is greater than or equal to the DC FIL DEGREE MAX LEVEL are tried to be brought up. We change the capacity levels of the correspondent open DCs provided that the move is not tabu. Regarding the computational results, we see that the best solution's improvement ratio after applying this move is around 0.0005.

4.2.9.
Optimizing the Open DCs' Orders Move. This move, which is called by optimizeDCsOrdersIf N ecessary sub-procedure is in fact a combination of moves and could only be applied when the intensification phase number is even. This move aims to optimize the orders and try to get free space in open DCs whose filling ratio is greater than or equal to the DC FIL DEGREE MAX LEVEL.

Diversification.
In the diversification (diversif y) procedure, we change the input solution's location-allocation decision until we get a new solution with different location-allocation decision. We check the triedLocAllocConf igsList to get a solution with different location-allocation configuration. When the new solution is created, this solution's location-allocation configuration is added to the triedLocAllocConf igsList.
To create new solutions, we apply one of the two moves to the input solution after comparing the randomly generated float number and the value of CLOSE OPEN DC RATIO parameter. In the first move which is applied when the randomly generated float number is less than or equal to the value of CLOSE OPEN DC RATIO parameter, we try to close one of the open DCs (I c ). While doing this operation we try to assign the retailers of this DC to the closed DC(s) that will be opened or to the other open DC(s). To make this decision, we compare a randomly generated float number and the value of the CLOSED DCS PRI RATIO parameter. If the generated float number is less than or equal to the value of CLOSED DCS PRI RATIO parameter, we try to assign the retailers to the new opened DC(s); otherwise, we try to assign them to the already opened DC(s). In the second move which is applied when the randomly generated float number is greater than the value of the CLOSE OPEN DC RATIO parameter, we select one of the opened DC(s) (I c ) and try to assign a fraction of its retailers to the other already opened DC(s) or close DC(s) which will be opened. The number of retailers that will be randomly selected from DC I c is found by multiplying the number of total assigned retailers and the ASSIGN RETA RATIO parameter. In these moves, before assigning the retailers to the already opened DCs or new opened DCs, we maximize the capacity levels of the DCs to make the retailer assignments effectively.
In the diversification moves, we select the open DC I c randomly using one of the two roulette-wheel selection methods. In the first method, we use the capacity usage information of open DCs (lower capacity usage means higher probability to be chosen), whereas in the second method, we use the estimated total cost of the open DCs found by the f indLocAllocDecisions procedure (higher estimated total cost means higher probability to be chosen). While choosing the DC(s) in both moves to be opened, we randomly generate a float number and compare it with the value of CLOSED DCS RAND OPENING RATIO parameter. If the randomly generated float number is less than or equal to the value of the CLOSED DCS RAND OPENING RATIO parameter, the DC(s) that will be opened are selected randomly; otherwise they are ordered by ascending according to their estimated total cost and then selected consecutively.
At the end of the diversification procedure; we minimize the capacity levels of all open DCs, optimize the ordering decisions of the open DCs whose filling ratio are greater than or equal to the DC FIL DEGREE MAX LEVEL consecutively through the updateOrdersW ithOrderingM odel sub-procedure and refind the routes of the open DCs to which new retailers are assigned or from which some retailers are shifted.

Computational study.
We coded the mathematical model and heuristic methods in Java programming language and solved the mathematical models through IBM ILOG CPLEX 12.8.0 solver. We keep PRESOLVE and HEURISTIC options of the solver open. Except for the results in Table 19, we took all runs on a PC with Intel Xeon E5-2640 v4 (32 2.40 GHz cores) processor, 64 GB RAM and 64 bit Windows 10 Professional operating system. 5.1. Instances. Since there are no available benchmark instances for the MILRP, we generated random instances whose parameters are based on real life information and the Guerrero et al. [7] study. We name the instances as I − J − H − N − ID, where I is the number of potential DCs, J is the number of retailers, H is the number of periods, N is the number of products and ID is the identification number for the instance.
We generate the demand of retailer j for product n at period t through Normal distribution: d jtn ∼ Norm (µ jn , σ jn ), where µ jn ∈ [5,15] and σ jn ∈ [0, 5]. We locate all the nodes (Producers, DCs and retailers) randomly in an area of size 100 x 100. The length of the i − j arc in km, dist ij , is

NINT (.) function approximates its
parameter to the closest integer value. We generate the unit inventory holding cost of product n for a single period t ∈ H randomly in the interval [0.1, 0.4] for DCs i ∈ I and in the interval [0.5, 0.9] for retailers j ∈ J. Besides; we generate unit obsolescence cost ε randomly using interval [0.01, 0.02]. Total unit holding plus obsolescence cost of node k ∈ K in period t ∈ H for product n ∈ N, h ktn , is h kt(t+1)n + . The capacity of DC i with capacity level e, W e i , is NINT (DCF W Q e . DCV W Q . totD / |I|) where DCF W Q e is constant capacity quotient for capacity level e, DCV W Q is dynamic capacity quotient and totD (total demand) is Constant cost quotients for vehicle types (V F CQ) are {1.5, 2, 2.5} respectively. We generate constant cost of using type b vehicle at least once (V F BU C) randomly in the interval [2000,4000]. Ordering cost of product n by DC i, OC in , is equal to the sum of DC ordering cost of that product plus length of the arc between product n's producer and DC i. We generate DC i's ordering cost of product n randomly in the interval [50, 100].

Preliminary tests.
According to the preliminary tests in which solution time and gap (gap between relaxed model's cost and optimal cost) analysis are made, we see that the proposed two valid inequalities for the M1 model and only the second proposed valid inequality for the M2 model are effective. Therefore, we added these valid inequalities to the relevant models. Regarding the P EN ALT Y F ACT OR, T ABU T EN U RE, and ALF A constants used in the vehicle routing algorithm, we take the same low and high values used in the Pasha et al. [15] study to generate and assign the random values to these constants. On the other hand, to find the appropriate values for the constants that substantially affect the performance and solutions' quality, we use "Nearly Orthogonal Linear HyperCubes Design" proposed by Cioppa and Lucas [3]. It is stated in this study that pairwise correlations between factors become nearly zero by using their robust design. we entered the low and high values for each factor into the given excel sheet and generated 16 different parameter sets to use their design. We run the routing algorithm 15 times for each parameter set and assessed all the parameter sets by considering their routing costs' (best, worst, average, and standard deviation). As a result, we find the appropriate values presented in Table  8. We use the same methodology described in the previous paragraph to find the appropriate values for the constants INI TEMP, FRE TEMP, COOLING SPEED, NUM OF ITERS AT ALL TEMPS, NUM OF CONS TEMPS SOL NOT IMP and NUM OF PHASES. These constants are used in the simulated annealing part of the intensification process. We take the same low and high values used in Zhang et al. [19] for the first five parameters, whereas take 5-20 values for the NUM OF PHASES parameter. After making the relevant analysis, we find the appropriate values shown in Table 9. After solving many instances and assessing the results, we come up with the appropriate probabilities for the application of intensification moves shown in Table  10 and the appropriate values to be assigned to the important constants (presented in Table 11) used in the diversification process. We assign the rounded-up value of |I | / 5 + |J | / 20 to the totalNumOfDiversifications parameter.  Table 11. Values assigned to the important constants used in the diversification process The value assigned to the relevant parameter changes with the phase currently in. b The value assigned to the relevant parameter is generated randomly in this interval.

5.3.
Results. In Section 2.1, we mentioned that our model for MILRP with heterogeneous fleet (M1 model) is based on the ILRP model proposed by Guerrero et al. [7]. We also pointed out that we made some modifications to the based model to decrease the number of decision variables and make the model easily understandable. To compare our model (M1 model) with a new MILRP model, which is totally based on the Guerrero et al. [7] (M1 G model), five 3-5-3-2 MILRP instances are solved optimally through ILOG CPLEX solver. The results are given in Table 12. The first column gives the name of the instance. Columns 2-3 show the computation time in seconds elapsed until the optimal solution is found with the M1 G and M1 models. We see that all the computation times obtained with the M1 G model are higher than the computation times obtained with the M1 model and average computation time for the M1 G model is nearly two times the average computation time for M1 model. These findings prove the effectiveness of our modifications to the M1 G model.
After doing preliminary experiments, we see that only the very small-sized instances could be solved through ILOG CPLEX solver with a time limit of 2 hours. In comparing the solver and the heuristic methods, 20 MILRP instances are solved. We create the instances in a way that we can see both optimal solutions, easy versus difficult instances and different instances concerning the number of potential DCs, retailers, periods and products. The results are given in Table 13. The first column gives the name of the instance. Columns 2-3 show the gap between lower and upper bound when instances are solved through the solver preset within 1200 s and 7200 s (2 hours). In column four, we see the cost of the best solution found within 7200 s. Column five gives the computation time in seconds elapsed until the best solution is encountered. Columns 6-8 present the average results for SHM in five runs. In these columns, we see the average cost of the solution, the average gap between the cost found by SHM and the cost found by the solver, and the average computation time in seconds. Columns 9-12 present the average results for HHM in five runs. In these columns, we see the average cost of the solution, the average gap between the cost found by HHM and cost found by the solver, the average gap between the cost found by HHM and the cost found by SHM and the average computation time in seconds.
Looking at the values in Table 13, we notice that only 5 out of 20 very smallsized instances could be solved exactly through the solver within a time limit of 7200 s. The quality of the solutions found with the solver gets worse quickly as the problem's size increases. This result proves the difficulty of the MILRP with heterogeneous fleet.
When we look at the results of the first ten instances, which are solved exactly or almost exactly, we see that average gaps between the cost found by the solver and costs found by the SHM and HHM are 7.4% and 6.1%, respectively. We can say that these gaps are acceptable. Besides, on average, the solutions found with SHM are 51.8%, and solutions found with HHM are 53.9% better than the solutions found with the solver regarding all instances. On average HHM outperforms SHM, but SHM outperforms the HHM prominently in two instances. We know that the reason behind this result is the exact inventory management decisions in SHM.
The average computing time for SHM is only 2.8 s, nearly one fifth of the average computing time for HHM (12.8 s). After making an in-depth analysis to understand the reason behind this significant difference, we noticed that SHM uses 99% CPU power of the PC (all (32) cores of the CPU) while solving the inventory management model through a multi-threaded CPLEX solver, which constitutes nearly whole part of the SHM computing time. On the other hand, HHM only uses 5% CPU power of the PC (only one core of the CPU). Nevertheless, computing times for both methods are acceptable.
The analysis regarding the solution time and solution quality of heuristic methods when solving the randomly generated small, medium and large-sized instances is shown in Table 14. The first column gives the name of the instance. Columns 2-5 present the average results of SHM run in five times. In these columns, we see the average cost, average computation time required to solve the instance, the average computation time needed to find inventory management decisions of the instance (IM CPU) and the average percentage of this time (IM CPU) in total computation time. Columns 6-8 present the average results for HHM in five runs. In these columns, we see the average cost of the solution, the average gap between the cost found by HHM and cost found by SHM and the average computation time in seconds.
When we look at the results for small-sized instances, we see that average computation time for SHM is 0.4 s and 73.4% (0.3 s) amount of this time is used to find the inventory management decisions. On the other hand, the proposed HHM spends 66.2 s for these instances and outperforms SHM by 4.4% on average. Considering the medium-sized instances, we notice that the HHM spends 1.5 times more computation time than SHM, but at the expense of this time, it outperforms SHM by 3.9% on average.
Finally, when we look at the results for large-sized instances, we see that the average computation time for SHM is 3572.6 s (nearly one hour) and the percentage of this time spent on inventory management decisions is similar to the mediumsized instances. On the other hand, HHM spends 1218.9 s (nearly 20 minutes) and outperforms SHM by 4.5% on average. Consequently, we can also say that these CPU times are acceptable for both methods considering the studied problem's difficulty.
After looking at the GAP values in Table 14, we notice that the most significant average gap between the compared solutions' costs occurs at 5-50-3-3-P5 instance. When we seek the reason behind this considerable difference and see that the best solutions found by SHM runs includes two open DCs, whereas the best solutions found by HHM runs includes only one open DC.
The analyses regarding the cost weight factors α, β, θ and γ used in the M1 model's objective function are presented in the following tables. We made all the analyses on medium-sized instances, which are solved five times through the HHM. Tables present the relevant values of the best solutions found. When we look at the results in Table 17, we see that when the value of the inventory cost weight factor for DCs changes to 100 times its original value, the total amount of products in open DCs' inventories decreases as expected. Besides, total number of orders made by open DCs increases to decrease the total amount of products in open DCs' inventories. TotInv OfDCs c total amount of products in retailers' inventories, the total number of distributions made to the retailers and the total amount of products in open DCs' inventories increases. We also notice that as a balancing effect, total number of orders made by open DCs decreases a little bit. Remember, in Section 3.3; we used SEH (shrinking and expanding) subprocedures in the f indRoutes procedure to shrink or expand the routes taken as input for diversification purpose. In Fig. 3 we see the effect of the SEH sub-procedures. We run the f indRoutes procedure on a randomly generated 5-50-3-3 instance and graph the found best solutions with and without the SEH sub-procedure.
As Fig. 3 shows, the f indRoutes algorithm starts with a good initial solution of 86590.8 but improves this solution slightly without the SEH sub-procedure. On the other hand, the f indRoutes algorithm starts with a relatively bad initial solution of 90762.8 and improves this solution many times with SEH sub-procedure's help. To see the effect of fixed fleet strategies used in the f indAllRoutes procedure, we randomly create five new 5-50-3-3 instances and run the SHM on these instances five times with and without fixed fleet strategies. The results are seen in Table 19. The second column show the best, the third column show the average distribution plus using vehicles at least once cost and the fourth column show the average computation time of the runs when the instances are solved without fixed fleet strategies. On the other hand, columns 5-7 show the same values when the instances are solved with fixed fleet strategies. We see that on the average, the average cost of the active case is nearly 13.5% lower than the passive case with the expense of almost two times of computation time. Computation times are short so this drawback is acceptable.
We present the objective function value of the current and global best solutions for one of the runs of the 5-15-3-3-P1 instance in Fig. 4 to validate the intensification and diversification procedures of HHM. In this run, the intensification procedure is called three times and all of the phases (12 phases) are recorded. On the other hand, the diversification procedure is called two times and records are seen above 13. and 26. x axis labels. 6. Conclusions. We present the MILRP with heterogeneous fleet for the first time as a comprehensive/generic approach to optimize many supply chain design scenarios including medicine/food delivery in humanitarian logistics, ammunition and other supplies distribution in the military logistics and cases in which leasing distribution centers are common, distribution to customers are made in different frequencies so the routing decisions change per period. MILRP with heterogeneous fleet considers a supply chain, which consists of multiple producers, potential distribution centers (DCs) with opening capacity levels, geographically scattered retailers, each of which has deterministic demand over a discrete planning horizon and heterogeneous fleet of vehicles. The goal is determining a set of DCs with their capacity levels to open (location decisions), assigning retailers to the opened DCs (allocation decisions) and for all periods finding the product quantities to be ordered by and distributed from DCs (inventory management decisions) and determining the fleet and routes (routing decisions) simultaneously to satisfy the demands of retailers with minimum cost.
We propose a mixed-integer linear programming model to describe the MILRP with heterogeneous fleet and strengthened it by adapting two valid inequalities proposed for ILRP. We see that only the very small-sized instances can be solved exactly within a reasonable computation time via ILOG CPLEX commercial solver after doing preliminary experiments. Therefore, we develop two heuristic methods (sequential and hybrid).
In the sequential heuristic method (SHM); we find the location-allocation decisions through a greedy algorithm first. Then, we solve our proposed inventory management model to find inventory management decisions exactly. Finally, we determine the vehicle routing decisions through a method based on tabu search metaheuristic. On the other hand, the hybrid heuristic method (HHM) consists of initialization, intensification and diversification steps. In the initialization step, we create the initial solution. Then, through the intensification process, which is the hybridization of simulated annealing and tabu search metaheuristics we improve this solution. Finally, with the diversification process, we change the locationallocation decisions of the best solution found by the intensification process and call the intensification process with this new solution. This loop is called for a pre-determined time and the best global solution is found. In order to create different solutions (having different location and allocation decisions), the location and allocation configurations tried before are stored in a list, which is updated by intensification and diversification processes.
Since there are no benchmark instances in the literature for MILRP with heterogeneous fleet, impact of two valid inequalities on both models (M1 model for MILRP and M2 model for inventory management decisions) and performance/solution quality of the developed heuristics and the ILOG CPLEX commercial solver are tested on randomly generated realistic instances. Results show that the two valid inequalities proposed for the M1 model and only the second valid inequality proposed for the M2 model are effective. In terms of solution quality, we can say that the quality of solutions found by SHM and HHM are acceptable for very small sized instances but we cannot guarantee this for the other sized instances due to the inability of making a comparison with solver. Nevertheless, we see that both heuristic methods provide important savings in acceptable run times compared to the commercial solver.
When we compare the heuristic methods in terms of solution time, we see that average computation time for SHM is substantially lower than HHM for small and medium sized instances, but it is the other way round for the large-sized instances. When we compare them in terms of solution quality, we see that HHM outperforms SHM in all sized instances.
The analyses made on medium-sized instances regarding the ordering, distribution and inventory cost weight factors, which are used in the M1 model show that all the results are compatible with the expected results. The presented benchmarks at the end of the results section prove some of our developed algorithms' effectiveness.
In conclusion, we can use both heuristic methods as a "what if" analysis tool for MILRP with heterogeneous fleet, which is a generic and frequently came across supply chain scenario. As future research, we are planning to improve our routing methodology first. We will then study the MILRP with two routing levels (routing added between producers and DCs) or the MILRP with time window constraints on distributions to retailers.