Solving a Novel Inventory Location Model with Stochastic Constraints and (𝑅,𝑠,𝑆) Inventory Control Policy

Abstract

We solve a novel inventory-location model with a stochastic capacity constraint based on a periodic inventory control (ILM-PR) policy. The ILM-PR policy implies several changes with regard to other previous models proposed in the literature, which consider continuous review as their inventory policy. One of these changes is the inclusion of the undershoot concept, which has not been considered in previous ILM models in the literature. Based on our model, we are able to design a distribution network for a two-level supply chain, addressing both warehouse location and customer assignment decisions, whilst taking into consideration several aspects of inventory planning, in particular, evaluating the impact of the inventory control review period on the network configuration and system costs. Because the model is a very hard-to solve combinatorial nonlinear optimisation problem, we implemented two heuristics to solve it, namely, Tabu Search and Particle Swarm Optimisation. These approaches were tested over small instances in which they were able to find the optimal solution in just a few seconds. Because the model is a new one, a set of medium-size instances is provided that can be useful as a benchmark in future research. The heuristics showed a good convergence rate when applied to those instances. The results confirm that decision making over the inventory control policy has effects on the distribution network design.

1. Introduction

Distribution network design (DND) is one of the most important problems for companies that distribute products to their customers. The problem consists of selecting specific sites to install plants, warehouses, and distribution centres, assigning customers to serving and interconnecting facilities by flow assignment decisions. DND is typically solved as part of a sequential approach that simplifies associated tactical and operational issues. Hence, the decisions that have been omitted are tackled only after the DND problem has already been solved. This means that decisions related to the location of warehouses and allocation of customers are made without taking into account operational aspects such as transportation and inventory. This situation leads to suboptimal designs because operational decisions are restricted to the current network design.

This paper considers a two-level supply chain, in which a single plant serves a set of warehouses, which in turn serve a set of end customers or retailers. Unlike traditional approaches and in accordance with the recent inventory-location literature, in this paper, we incorporate the inventory control policy as a relevant factor that directly affects DND. A distinctive feature of our model is that it is based on a periodic inventory control policy (𝑅,𝑠,𝑆) for each distribution centre (DC) in a single-product scenario, which is important for those industries in which a continuous revision policy is not possible. Because our model is a very hard-to-solve nonlinear one, we have solved it using well-known heuristic approaches, namely, Tabu Search (TS) and Particle Swarm Optimisation (PSO). Moreover, as the model was recently developed by us, there is no set of instances that can be used as a benchmark. Thus, in this paper, we present two types of medium-sized instances. The first type considers uniformly distributed customers. The other one considers two clusters of customers. Warehouses are uniformly distributed in both types of instances.

The remainder of the paper is organised as follows. First, a literature review of the main ILM is presented in Section 2. Then, in Section 3, we present and analyse the stochastic capacity constraint under periodic review (𝑅,𝑠,𝑆). At the end of that section, the formulation of the model is presented. Section 4 briefly presents a review of both TS and PSO algorithms and some features of the corresponding implementations. Section 5 starts explaining the procedure used to generate the set of instances presented in this paper. In Section 5.1, the obtained results are presented. Finally, in the last section, some conclusions based on the numerical results are outlined.

2. Literature Review

Many authors have been focused on the DND problem and on ILM in particular during the last 15 years. For instance Simchi-Levi and Zhao [1], Mourits and Evers [2], Bradley and Arntzen [3], Miranda and Garrido [4], and Miranda [5] all analysed the levels of decision making related to DND and supply chain management (SCM). Daskin [6], Simchi-Levi et al. [7], and Drezner and Hamacher [8] present detailed FLP reviews and analysis. However, the traditional structure of FLP is not useful for considering interactions between facility location and inventory decisions or the impact of the latter on network configurations. For example, the risk pooling effect states that the total system safety stock is reduced when customers are served by a smaller number of warehouses. Daskin et al. [9] and Shen et al. [10] incorporate an (𝑄,RP) inventory control policy into the widely studied uncapacitated facility location problem (UFLP), establishing a safety stock at each site. Daskin et al. [9] employ Lagrangian relaxation to solve the model, whereas Shen et al. [10] reformulate the model as a set-covering problem and solve it using a column generation method. Based on the same inventory control policy, Miranda and Garrido [4] consider the order quantity for each warehouse as a decision variable and the capacitated facility location problem (CFLP) as a base framework. Finally, in Miranda and Garrido [11] and Ozsen et al. [12] authors handle capacity constraints by using previous inventory-location models.

More recently, Kumar and Tiwari [13] have presented an ILM that incorporates the risk pooling effect for both safety stock and running inventory. Additionally, in their model, the authors consider the effect produced when warehouses and end customers work jointly. Tancrez et al. [14] have presented a three-level supply chain nonlinear ILM that integrates location, allocation, and shipment sizes. In [15], the authors developed a model for the DND problem that considers the short lifetime of perishable products. To solve their model, the authors implemented a Lagrangian relaxation. One paper that has addressed the DND problem using periodic 3 inventory review is presented in [16]. In their paper, the authors consider a (𝑅,𝑆) inventory policy, which is slightly different from the one that is considered in this work.

The authors have presented different techniques to solve these models. For instance, Bard and Nananukul [17] proposed a branch and price algorithm for an integrated production and inventory routing problem. In [18], Lagrangian relaxation was used to solve a mixed integer linear programming model for multiple echelon and multiple commodity supply chain network design. Heuristics have also been used to solve DND problems. For instance, Armentano et al. [19] TS is used to minimise the production and inventory cost in a model that integrates production and distribution decisions by considering the capacity constraints of the plant. Askin et al. [20] implemented an evolutionary algorithm to solve an ILM considering multicommodity and distribution planning decisions. In their paper, the authors present a very comprehensive description of their genetic algorithm. To the best of our knowledge, PSO has not been considered to solve ILM problems.

As we stated previously, in this paper, we have considered a two-level supply chain with a single plant, a set of warehouses, and a set of end customers or retailers. Because we incorporate the inventory control policy as a relevant factor that directly affects DND in this paper, based on a periodic inventory control policy (𝑅,𝑠,𝑆) for each Distribution Center (DC), this model can be considered as a variant of the models presented previously in the literature [4, 9, 12, 21], which considers a policy of continuous inventory review (𝑄,𝑅). In next section, the new model proposed in this work is presented.

3. Inventory-Location Model with Stochastic Capacity Constraint

The model presented in [11] optimises warehouse location and customer assignment decisions, taking into account fixed installation, transportation, inventory, and fixed ordering costs. The authors assume that each warehouse 𝑖 operates a continuous inventory review policy based on a (𝑄,RP) policy to meet a stochastic demand, with mean 𝐷 (units of product per time unit) and variance 𝑉. It is also assumed that the plant takes a lead time, LT, to fill incoming orders from warehouse 𝑖. Stochastic constraint on inventory capacity is proposed, assuming a maximum inventory level for each warehouse 𝐼^cap. This constraint is based on chance constrained programming [22] and it fixes a maximum probability, 𝛽, of violating inventory capacity at peak times, which occurs only when orders arrive at the warehouses. This inventory level corresponds to the reorder point RP, which is stated in order to satisfy demand during LT with at least a probability of 1−𝛼, minus stochastic demand during LT, plus order quantity, RP−SD(LT)+𝑄. Thus, the inventory capacity constraint can be written as a deterministic nonlinear constraint as follows: 𝑄+𝑍_1−𝛼+𝑍_1−𝛽⋅√LT⋅√𝑉≤𝐼^cap. (1)

It may be noted that a similar and more conservative capacity constraint is proposed in [12], which ensures that inventory capacity is observed in 100% of the cases. However, when periodic review is considered, particularly assuming a (𝑅,𝑠,𝑆) inventory control policy [23], capacity constraint cannot be stated at any moment and does not take the same form as in (1). In a (𝑅,𝑠,𝑆) inventory control policy, inventory levels are reviewed only every 𝑅 period, and if the inventory level is lower than 𝑠, then an order is submitted to reach the objective level 𝑆. Consequently, order size must consider the well-known undershoot magnitude (US), which is the amount of items required in addition to 𝑠 to reach 𝑆 units of inventory. The average US as a function of demand mean and variance, 𝐷 and 𝑉, and for a review period 𝑅, can be computed as US()=𝑉𝐷,𝑉+2⋅𝐷𝑅⋅𝐷2. (2)

In terms of inventory capacity constraints, peak inventory levels are not controlled at any time, but only at specific times for each review period. The peak inventory level is reached only when orders arrive at the warehouse LT time units after the last order and naturally only if an order was submitted to the central warehouse or plant. Consequently, when an order arrives at a warehouse, the inventory level is (𝑠−US)_{Inventorylevelwhenanorderissubmitted}+(𝑆−𝑠+US)_{Submittedordersize}−SD(LT)_{StochasticdemandduringLT}=𝑆−SD(LT)_{Maximuminventorylevelwhenanorderarrives}. (3)

This expression is not surprising, as when an order is submitted to the plant, it is necessary that the total inventory position (on hand plus on order inventory) reaches level 𝑆 and that LT time units later, the inventory level is reduced by the demand during the LT, SD(LT). Similar to [11], in this paper, we propose that this constraint must be observed for each peak inventory instant (i.e., for each order period) with a fixed and known probability of 1−𝛽, but now assuming a periodic review as follows: 𝑆≤𝐼^cap+𝐷⋅LT−𝑍_1−𝛽⋅√𝑉⋅LT. (4) We then define the minimum order size, 𝑄, as 𝑆=𝑠+𝑄. (5) Consequently, constraint (4) can be rewritten as follows: 𝑄+𝑠≤𝐼^cap+𝐷⋅LT−𝑍_1−𝛽⋅√𝑉⋅LT. (6)

Finally, the reorder point 𝑠 is set to ensure that, for each time an order is not submitted (inventory level larger than 𝑠), the inventory level is large enough to fill the demand until the next order arrives in 𝑅+LT time units, with a probability or service level of 1−𝛼: 𝑠=𝐷⋅(LT+𝑅)+𝑍_1−𝛼⋅√𝑅+LT⋅√𝑉. (7)

Finally, substituting (7) into (6), the inventory capacity constraint is 𝑄+𝐷⋅𝑅+𝑍_1−𝛼⋅√𝑅+LT+𝑍_1−𝛽⋅√LT⋅√𝑉≤𝐼^cap. (8)

According to the previous inventory control assumption, we describe the proposed Inventory-Location Model with Stochastic Constraints on Inventory Capacity under Periodic Review (ILM-SCC-PR) as a stochastic nonlinear mixed-integer programming model (SNL-MIP). Based on a periodic (𝑅,𝑠,𝑆) inventory control policy, the safety stock to be included in the objective function is the average inventory level just before an order arrives at the warehouse: (𝑠−US)−𝐷⋅LT=𝐷⋅𝑅+𝑍_1−𝛼⋅√𝑅+LT⋅√𝑉𝑉−+2⋅𝐷𝐷⋅𝑅2. (9)

Additionally, inventory and ordering costs related to order size or inventory cycle are evaluated in terms of the minimum order size 𝑄 as a decision variable, as in the well-known EOQ model: OC⋅𝐷(𝑄+US)+HC⋅(𝑄+US)2. (10) The variables and parameters considered in the mathematical formulation are as follows.

Parameters𝑁: Number of available sites to install warehouses  𝑀: Number of customers that must be served by the installed warehouses  𝐹_𝑖: Fixed location cost of a warehouse on site 𝑖 ($/day)  𝐶_𝑖𝑗: Transportation unit cost from the warehouse on site 𝑖 to customer 𝑗 ($/unit)  OC_𝑖: Fixed ordering cost at site 𝑖 ($/order) HC_𝑖: Holding cost per time unit at site 𝑖 ($/day) 𝜇_𝑗: Mean demand of customer 𝑗 per day 𝜎_𝑗: Standard deviation of the demand of customer 𝑗 per day 𝐼^cap_𝑖: Capacity at warehouse 𝑖 𝑄^cap_𝑖: Maximum order size at warehouse 𝑖 𝑅_𝑖: Inventory check period at warehouse 𝑖 in days LT_𝑖: Average leadtime at warehouse 𝑖 in days 𝑍_1−𝛼: Value of the Standard Normal Distribution, which accumulates a probability of 1−𝛼 𝑍_1−𝛽: Value of the Standard Normal Distribution, which accumulates a probability of 1−𝛽.

Variables  𝑋_𝑖: Binary variable. It is equal to 1 if a warehouse is installed on site 𝑖 and 0 otherwise 𝑌_𝑖𝑗: Binary variable. It is equal to 1 if warehouse i serves customer 𝑗 and 0 otherwise 𝑄_𝑖: Order size at warehouse 𝑖. It is greater than 0 if 𝐷_𝑖>0, and 0 otherwise 𝐷_𝑖: Total mean demand at warehouse 𝑖. It is greater than 0 if there exist at least one 𝑌_𝑖𝑗>0 and 0 otherwise 𝑉_𝑖: Total variance of the demand at warehouse 𝑖. It is greater than 0 if 𝐷_𝑖>0 and 0 otherwise.Then, the optimisation problem is Min ^𝑁_𝑖=1𝐹_𝑖⋅𝑋_𝑖+^𝑁𝑀_𝑖=1𝑗=1𝐶_𝑖𝑗⋅𝑌_𝑖𝑗+^𝑁_𝑖=1OC_𝑖⋅𝐷_𝑖𝑄_𝑖+US_𝑖+HC_𝑖⋅𝑄_𝑖+US_𝑖2+^𝑁_𝑖=1HC_𝑖𝐷_𝑖⋅𝑅_𝑖+𝑍_1−𝛼⋅𝑅_𝑖+LT_𝑖⋅𝑉_𝑖−US_𝑖 (11) s.t.: ^𝑁_𝑖=1𝑌_𝑖𝑗=1∀𝑗=1,2,…,𝑀,𝑌_𝑖𝑗≤𝑋_𝑖∀𝑖=1,2,…,𝑁,∀𝑗=1,2,…,𝑀,𝑄_𝑖+𝐷_𝑖⋅𝑅_𝑖+𝑍_1−𝛼⋅𝑅_𝑖+LT_𝑖+𝑍_1−𝛽⋅LT_𝑖⋅𝑉_𝑖≤𝐼^cap_𝑖⋅𝑋_𝑖∀𝑖=1,2,…,𝑁,𝑄_𝑖+US_𝑖≤𝑄^cap_𝑖⋅𝑋_𝑖∀𝑖=1,2,…,𝑁,𝐷_𝑖=^𝑀_𝑗=1𝜇⋅𝑌_𝑖𝑗∀𝑖=1,2,…,𝑁,𝑉_𝑖=^𝑀_𝑗=1𝜎²_𝑗⋅𝑌_𝑖𝑗∀𝑖=1,2,…,𝑁,US_𝑖=𝑉_𝑖2⋅𝐷_𝑖+𝑅_𝑖⋅𝐷_𝑖2if𝐷_𝑖0>0otherwise𝑋_𝑖,𝑌_𝑖𝑗∈{0,1}∀𝑖=1,…,𝑁;∀𝑗=1,…,𝑀. (12)(13)(14)(15)(16)(17)(18)(19)

Equation (11) is the total system cost. The first term is the fixed setup and operating cost when opening warehouses. The second term is the daily transport cost between warehouse and customers. The third term contains fixed order and inventory costs related to warehouse order size. The fourth term represents the storage cost associated with safety stock at each warehouse. Equation (12) ensures that retailers are served by only one warehouse. Inequality (13) ensures that customers are only assigned to installed warehouses (𝑋_𝑖=1). Inequality (14) ensures that the inventory capacity of each warehouse is respected at least with a probability 1−𝛼. Inequality (15) ensures that the order size is below the maximum order size 𝑄^cap_𝑖 allowed to warehouse 𝑖. Equations (16) and (17) compute the mean and variance of the total demand served by each warehouse 𝑖, respectively. Equations (18) determine the value of US_𝑖. Finally, (19) states integrality (0−1) for the variables 𝑌_𝑖𝑗 and 𝑋_𝑖. This model is NP hard because it is clearly an extension of the UFLP, which is already NP hard (UFLP can be obtained just by using HC=0 and OC=0). In addition, the objective function and one constraint are nonlinear, resulting in a model that is very difficult to solve to optimality using classical mathematical programming techniques. Therefore, in this paper, we attempt to solve medium-size instances by mean of two well-known heuristics called Tabu Search (TS) and Particle Swarm Optimisation (PSO).

4. Heuristic Methods

Heuristic methods are a common approach to solve hard combinatorial optimisation problems. Despite the fact that heuristics do not guarantee optimality, the solutions provided by them can be considered as good suboptimal ones. In contrast exact methods guarantee optimality; however, they usually fail when dealing with medium- and large-sized problems. In this paper, two heuristics have been separately considered to solve our DND problem. The first one corresponds to a well-known local search called Tabu Search. The second one is an evolutionary algorithm called Particle Swarm Optimisation. In the next subsections, an overview of these two heuristics is provided.

4.1. Tabu Search

We can describe the TS approach as a local search technique guided by the use of adaptive or flexible memory structures. However, such a general definition fails to show the specificity of TS and could be confused with other types of Greedy Random Adaptive Search Procedure (GRASP) algorithms. The variety of the tools and search principles introduced and described in [24] are such that the TS can be considered as the seed of a general framework for modern heuristic search [25]. TS has been applied to several combinatorial optimisation problems (see, for instance, [26–29]). So too other techniques [30–32] have been used to solve different variants of the FLP. The TS is essentially a local search algorithm; that is, it needs to “exchange information” with its neighbours. To do that, first, the neighbourhood must be defined. Typically, TS algorithm has only one neighbourhood move that defines a set of possible candidate solutions. This move is used across all executions of the TS algorithm. Because different neighbourhood definitions can lead to different results in this paper, we have implemented two types of neighbourhoods. The first one is defined by a change in the allocation of a customer from one facility to another without any restriction; that is, the new facility could be either open or closed at the allocation moment. The second one is defined by a change in the allocation of a customer from one facility to another in which the new facility was opened previous to the movement 90% of the time; in other words, the probability of moving a customer to a closed facility is equal to 10%. The TS heuristic provides a diversification mechanism that allows it to get out of low-quality neighbourhoods and “jump” to explore new neighbourhoods. The diversification mechanism implemented in this study is a restart method, which re-initialises a current solution without losing the best solution found by the algorithm. As stated above, the restart method is used in case the TS is unable to move out of low-quality neighbourhoods. The TS algorithm requires the following structures to be implemented. (i)solutionVector: a vector with a size equal to the number of customers. For each customer this vector contains the DC to which the customer is allocated. (ii)bestSolutionVector: a vector corresponding to the best solutionVector found during the TS execution. (iii)candidateSolutionVector: a vector resulting from application of the neighbourhood movement to a solutionVector. (iv)candidateList: a list of candidateSolutionVectors. From this list, we obtain the next solutionVector, which corresponds to the best candidateSolution from the list. (v)tabuList: an 𝑁×𝑀 matrix that contains those neighbourhood movements that are prohibited in the current TS iteration. Additionally, our TS implementation also requires the following parameter list. (i)iterationNumber: Total iteration number to end the algorithm. (ii)diversificationBound: total iteration number to apply diversification criterion (restart method). (iii)tabuListSize: size of tabuList. This means the number of iterations for which a specific movement remains banned. (iv)listOfCandidates: size of the set of candidateList.

The TS algorithm starts with a random solution in which each customer is allocated to a randomly selected DC. When all customers are assigned, TS validates the model constraints. If the constraints are satisfied, the initial solution (𝑆₀) is assigned as a best solution (𝑆_best). A list of candidate solutions (𝑁𝑆₀) is generated after that. They are neighbours of 𝑆₀. The size of the 𝑁𝑆₀ set corresponds to the listOfCandidates parameter. Once the set of 𝑁𝑆₀ has been generated, the TS algorithm selects the best candidate solution (𝑁𝑆^best₀), which should not stay on the tabu list. If the 𝑁𝑆^best₀ actually is in the tabuList, it cannot be selected unless its cost is less than the cost of the 𝑆_best found until this moment (aspiration criterion). When a new current solution 𝑆₁ is selected, the tabuList is updated, and a new iteration starts. The algorithm stops when the number of iterations 𝑘 is equal to the parameter iterationNumber. Below, we can see the general framework of our TS implementation.

4.2. Particle Swarm Optimisation

The PSO algorithm was proposed and developed by Kennedy and Eberhart [33–37]. Although the PSO algorithm was initially designed to tackle continuous optimisation problems, during the last decade, it has proven to be a very good alternative for solving combinatorial optimisation problems. In fact, several articles have used this technique to solve complex combinatorial optimisation problems; see, for example, [38, 39]. For a comprehensive analysis of publications concerning PSO approaches, see [40].

Unlike other evolutionary algorithms, PSO does not use the “survival” concept. This is because all particles are kept “alive” throughout the algorithm execution time and at no time is their survival threatened. The algorithm starts when a set of particles 𝒮 (or swarm) is initialised. Each particle 𝑝^𝑘∈𝒮, with 𝑘=1,…,𝑃, is represented by a vector in which the value of each element 𝑝^𝑘_𝑗 with 𝑗=1,…,𝑀 corresponds to the warehouse that has been set to a customer 𝑗. In our algorithm, the initialisation is a random process. For each particle, we know its current position (denoted by 𝑝^𝑘) and its previous best position (denoted by 𝑝^𝑘(best)). We also calculate for each particle the best current neighbour (denoted by 𝑝^{𝑘(neighbour)}) and the best particle (denoted by 𝑝^best) found so far. The position of each individual particle is adjusted at each iteration by considering its previous best positions, the position of the best neighbour, and the position of the global best solution found so far (see Algorithm 1).

begin    𝑘 = 0;   Generate initialSolution 𝑆𝑘;    𝑆best = 𝑆𝑘;    while 𝑘<  iterationNumber do    NS𝑘 = Generate Neighbourhood (𝑆𝑘);    find best (NS𝑘);    while    NSbest𝑘  isTabu   do       if   cost of NSbest𝑘 < cost of 𝑆best   then         aspirationCriterion;      reeplace NSbest𝑘;    If   cost of nsbest𝑘 < cost of 𝑆best   then       𝑆best=nsbest𝑘;       iteration without improvements = 0;    iteration without improvements + 1;    𝑆𝑘=nsbest𝑘;    update (tabuList);    check (diversificationCriterion);    𝑘=𝑘+1 end

Algorithm 1 

General algorithmic frame for tabu search.

The general frame for our PSO algorithm is presented below.

The most important step in Algorithm 2 is the particle updating. As we have said before, the PSO algorithm was designed for continuous optimisation problems. Therefore, we need to adapt some of its steps to solve combinatorial problems. One important change is in the strategy used to move particle 𝑝^𝑘 from its position in the iteration 𝑡 to its next position in iteration 𝑡+1. To do that, we have considered three possible alternatives: the first one (namely 𝑠₁) corresponds to a change in the allocation of one of the customers (randomly selected) to a new (randomly selected) warehouse. The second alternative (𝑠₂) is a change in the allocation of one of the customers (randomly selected) to a new warehouse following the allocation of either the bestNeighbour (𝑝^{𝑘(neighbour)}) or the prior best position (𝑝^𝑘(best)). Finally, the third alternative (𝑠₃) corresponds to the same movement but now considering the allocation from the global best solution reached so far 𝑝^best. Therefore, the position of particle 𝑝^𝑘 in the iteration 𝑡+1 will be the best particle among 𝑠₁, 𝑠₂, and 𝑠₃. Below we present the pseudocode for our particle updating process (see Algorithm 3).

begin    Initialise set 𝑆;    Calculate 𝑝best;    while   maxIteration  not  reached   do      foreach    𝑝𝑘∈𝑆   do          Calculate 𝑝𝑘(best);          Calculate 𝑝𝑘(𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑢𝑟);      Update 𝑝𝑘;      Calculate 𝑝best; end

Algorithm 2 

General algorithmic frame for particle swarm optimisation.

begin    Particles 𝑠1, 𝑠2, 𝑠3;    If   inertia   >   𝑟𝑎𝑛𝑑[0,1]   then      𝑠1= LocalMove𝑝𝑘;    If   cognitiveFactor   >   𝑟𝑎𝑛𝑑[0,1]   then      𝑠2= LocalMove𝑝𝑘(best);    𝑠2= LocalMove𝑝𝑘(𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑢𝑟);    If   globalCognitiveFactor   >   𝑟𝑎𝑛𝑑[0,1]   then      𝑠3= LocalMove𝑝best;    𝑝𝑘= best𝑠1,𝑠2,𝑠3; end

Algorithm 3 

Particle update.

As with local search algorithms, the PSO approach can also get trapped in local optima if the global best and local best positions are equal to the position of the particle over a number of iterations [41]. One distinctive feature of our algorithm is the neighbourhood definition. The idea is that we know that every particle 𝑝^𝑘 is able to see the best particle at the current iteration. Then, we assume that each particle 𝑘 can see, at least, any other particle 𝑙 for which the distance 𝑑_𝑘−𝑙 from 𝑝^𝑘 to 𝑝^𝑙 is less than or equal to distance 𝑑_𝑘-best from 𝑝^𝑘 to 𝑝^best. Distance 𝑑_𝑘−𝑙 is calculated as follows: dist_𝑘−𝑙=^𝑀_𝑗=0𝑝^𝑘_𝑗−𝑝^𝑙_𝑗². (20)

It is clear that in this paper, we do not attempt to develop an original strategy to solve our problem. Instead, we are proposing very simple procedures to obtain solutions for our new model that can be used as a baseline in future research. Moreover, validating some of the main model assumptions is also one of the goals of this work.

5. Experiments and Computational Results

In this section, we present the experiments we have carried out, and we draft some conclusions about our numerical results. As we have noted before in this papaer, our goal is not to state whether the TS or PSO approach is the best way to solve this problem. In fact, we firmly believe that other approaches might be more effective than our approach, and, therefore, applying different strategies to solve our model appears to be a fruitful area for future research into artificial intelligence field. In spite of this, due to the simplicity that both TS and PSO offer, we present these results as a baseline for future studies. We have selected these two algorithms as samples of local search and evolutionary strategies, respectively. Therefore, two of the most important approaches used to solve large and complex combinatorial optimisation problems are covered with these two algorithms.

Moreover, in this section, we present a set of instances that can be used as a benchmark in future studies. The benchmark set consists of two base instances, namely, 𝐶₁ and 𝐶₂, which means that customers are distributed uniformly and in clusters, respectively. Therefore, for instances of type 𝐶₁ customers and DCs are randomly distributed, whereas, for instances of type 𝐶₂, only DCs are distributed over the entire area. Customers are distributed randomly along a specific perimeter around the centre of the cluster. The maximum distance allowed between a customer and the centre of the cluster is a parameter. In our case, the number of clusters was equal to 2 and the ratio was 25[km]. These two base instances (𝐶₁ and 𝐶₂) were used for all the corresponding instances. Customer demand 𝑑_𝑗 is calculated using a uniform distribution 𝑈[30,100]. Transportation costs (TC) 𝑌_𝑖𝑗 correspond to 1/5 of the Euclidean distance between a customer 𝑗 and a warehouse 𝑖. Fixed costs (FC) are calculated as the sum of a constant value 𝑏 (in our case 5000) and a factor dependent on the distance between the warehouse and its closest cluster (10×𝑏/mindist). With this expression we make those facilities that are close to the “city centre” more expensive. In the same way, the capacity of the warehouse 𝐼^cap will decrease as the distance to the centre of the cluster becomes small. 𝐼^cap_𝑖 is calculated as the sum of a base value 𝑐=100 plus 10×mindist. The order cost (OC) and holding cost (HC) depend on the 𝐼^cap. The larger the capacity, the smaller the inventory costs. LT is a random number between 2 and 4. The size of the order 𝑄^cap is determined as the sum of a constant equal to 100 plus a percentage of the capacity of the warehouse 𝐼^cap. In our case, this percentage was 50%. FC, HC, OC, Service Levels (𝑍_1−𝛼), LT, 𝐼^cap, and 𝑄^cap for each DC are the same for the two base instances. Figure 1 shows these configurations.

(a)

(b)

Figure 1 

Customer and warehouse distributions. (a) Shows that customers (blue squares) are uniformly distributed. (b) Shows that customers are distributed over two clusters (blue squares and green triangles). In both cases, warehouses (orange diamonds are uniformly distributed).

From these two main configurations, we generated a set of 22 instances based on them (11 for each one). To do that, we modified values of FC, HC, and OC in ±25%, and additionally we have varied parameter 𝑅 among values 1, 2 and 3. Then, each instance attempts to evaluate the impact on the network configuration (and on its costs) produced by either increasing or decreasing the parameter values, one at a time. As base instances, we have 𝐶₁ and 𝐶₂ and based on them we have a set of instances such as 𝐶₁FC_.75, which means that the FCs from base instance 𝐶₁ have been reduced to 75% of their original value.

5.1. Computational Results and Discussion

The computational experiments were performed on an Intel Core Duo processor CPU T2700 2.33 GHz with 2 GB of RAM, using the Ubuntu 12.04 operating system. Both the TS and PSO algorithms were implemented in the JAVA programming language using NetBeans IDE. To validate the algorithms and measure their convergence, we used a set of small instances from [42] that have an optimal solution that was obtained by means of an enumerative method. This enumerative algorithm was executed for approximately 4 hours for each instance, whereas the heuristic algorithms took less than 5 seconds to solve each one. Furthermore, they found the optimal solution all 10 times that they were executed over each instance. We also ran a random search algorithm to determine the real improvement provided by our heuristic. On average, savings when using heuristics algorithms were approximately 300% in all instances.

We fixed both the nonstockout probability and the nonviolate capacity constraint probability at 95%. Hence, the values for 𝛼 and 𝛽 are fixed at 1.648. As we mentioned in Section 4 of this paper, we implemented two types of neighbourhood movements for our TS algorithm. Table 1 shows a summary of the results obtained using both movements for all 22 instances.

Table 1 

Summary results. Comparison between the results obtained by movements 𝑚₁ and 𝑚₂.

Column Instance corresponds to the evaluated instance. Columns 𝑥_𝑚1 and 𝑥_𝑚2 correspond to the average value obtained after executing the TS algorithms 10 times using movements 𝑀1 and 𝑀2, respectively. In the same way, columns 𝜎_𝑚1 and 𝜎_𝑚2 correspond to the standard deviation of the results. Finally, columns Time_𝑚1 (sec) and Time_𝑚2 (sec) show the average time needed by our TS implementation to find the best solution.

Table 1 shows how, predictably, the average cost increases with the value of period review 𝑅. Furthermore, the DND is also affected by the change in the value of 𝑅. In several instances, when 𝑅 increases, the number of open DCs follows it. Figures 2(a) and 2(b) illustrate this situation.

(a)

(b)

Figure 2 

(a) Shows the DND obtained when parameter 𝑅=1 and (b) shows the DND obtained when parameter 𝑅=3.

These two situations are produced for the same reason: when the value of 𝑅 increases, the safety stock must be greater to avoid stockout during the period without review. Thus, DCs cannot attend the same number of customers when the safety stock is increased because of their limited capacity. For this reason, the model is forced to open an extra DC to satisfy total customer demand. This relation is especially important because it shows that the value of 𝑅 not only affects the cost of the distribution network but also its design.

We now compare the best solutions obtained by the TS algorithm with the ones obtained by the PSO method. Table 2 shows the results of this comparison.

Table 2 

Comparison between results obtained by the local search (TS) and the evolutionary algorithm (PSO).

According to the results the PSO algorithm seems to perform slightly better than the TS. That is more evident when we look at 𝐶₂ instances where PSO shows a better performance in 9 out of 11 cases. On the other hand, TS is clearly better than PSO when 𝐶₁ instances are considered. These are quite interesting results because they suggest that we should take into account how customers are distributed to decide among different heuristic techniques. It is clear that further research is needed in this field.

Because we do not have any bound for our instances, we cannot say anything about the quality of our solutions. Another interesting research line would be to find lower bounds for this model to determine how good the solutions provided by the heuristic methods are.

Moreover, when we examine the components of our objective function, we can see how they behave. In Table 3, columns FC (%), TC (%), Inv (%), and SS (%) show the change (on average) of Location, Allocation, Inventory, and Safety Stock cost as a percentage of the base instance 𝐶_𝑖𝑅₁.

Table 3 

Variation of cost distribution at each instance as a percentage of base instances 𝐶_𝑖𝑅₁.

Another interesting finding is that individual subcosts (FC, TC, Inventory, and SS) do not increase to the same extent that parameters do. In other words, when, for example, FC is increased (decreased) by 25%, the portion of the total cost corresponding to FC does not increase by 25% when compared with its corresponding base instance. This can be associated with the ability of the model to seek other location/allocation alternatives to keep the total cost as low as possible.

Regarding our algorithms, they showed good performance in terms of time and convergence. As we mentioned before, we cannot compare our results with other works because this is a very new model. Figure 3 shows the convergence of the TS algorithm for a specific instance. Here, we can see how the diversification method was invoked several times (peaks in the solid line) during the execution after finding a (presumably) local optima.

Figure 3 

Current solution costs. The peaks correspond to the values of the current solution when the restart method is invoked.

Despite the fact that restart method proved to be effective in providing some diversification to our algorithm (and consequently allowing us to move away from local optima), other diversification strategies can be implemented to exploit information from previous iterations and, consequently, provide the algorithm a higher exploration level. With regard to the aspiration criterion, this was invoked 30 times on average for each execution with a lower bound of 4 and an upper bound of 71. These values show how the frequency-based memory of the TS works.

6. Conclusions and Future Work

In this paper, we solve a novel inventory-location model that integrates inventory decisions at the strategic level. As a first contribution, our model considers a periodic inventory review policy (𝑅,𝑠,𝑆) unlike previous models that assume continuous inventory review policies. We have also generated a set of 22 medium-sized instances for both uniformly and 2-clustered distributed customers. These medium-sized instances consider 50 DCs and 100 customers. This set of instances was solved using both TS and PSO heuristics independently which yielded an approximation to the optimal solution for each instance. Despite the simplicity of both approaches, good solutions (w.r.t. randomly generated ones) were found in an acceptable amount of time. The results obtained confirm the importance of the value of 𝑅, which affects both the cost and configuration of the distribution network. A sensitivity analysis over the key parameters of the model was performed. The numerical results showed the effect on the network produced by changes in both inventory cost (HC and OC) and TC. As future work it is possible to integrate other tactical elements such as routing or layout decisions to our ILM, to make them even more representative of reality and therefore improve the quality of the solutions. Examples of these elements are the transportation decisions or the DC layout decisions. Based on the obtained results, we can state (preliminary) that evolutionary algorithms and particularly those based on swarm intelligence such as PSO should be implemented to obtain better solutions for those instances that involve clusters. Local search strategies should be considered for those instances with uniformly distributed customers. However, further investigation is needed in this area.

Moreover, some matheuristics approaches (ones that mix mathematical programming and metaheuristics) seem to be another interesting research area to explore. Finally, it would also be interesting to calculate some lower bounds for this model in order to measure how far from the optimal solution our approximations are.

Acknowledgment

The authors are indebted to anonymous referee comments and the editor for their valuable comments through the review process. The paper remarkably improved through their recommendations; yet, the authors are responsible for any remaining errors. Fernando Paredes is supported by FONDECYT-Chile Grant 1130455.