1. INTRODUCTION
A supply chain is a network of activities associated with the flow of goods, starting from the conversion of raw materials into products and ending with the delivery of final products to consumers at the right quantity, time, and place at the lowest possible cost (maximum profit) and with a satisfactory service level (Pahlevan et al., 2021;Amiri, 2006;Javid and Azad, 2010). Supply chains often comprise a series of facilities including manufacturing plants, depots, distribution centers, service centers, and retailers (Bektas, 2006;PérezRodríguez, 2021).
One of the key questions in the design of logistics systems is where to establish distribution centers. Manufacturers can greatly reduce their production and distribution costs by finding the best locations for their distribution warehouses and the best routes for their transportation operations (BenTal et al., 2009;Guillén et al., 2005). These two objectives can be pursued simultaneously by combining facility location and vehicle routing problems. In this paper, this approach is used to optimize the routes by which vehicles transport products simultaneously with the location of warehouses whereby product distribution is conducted; an approach that is expected to provide better solutions than optimizing each aspect independently.
The limited shelf life of perishable goods imposes new challenges on supply chain management. Obviously, this is because perishable goods that remain in inventory longer than their shelf life need to be separated and discarded. Therefore, the quantities sent to retailers should be limited based on shelf life as well as the capacity of retailers.
In a study by Qian et al. (2020), they developed a model for singlesource plant location. Harris et al. (2020) (2021) also studied the subject of plant location for a given set of customers. Pirkul and Jayaraman (1998) provided a model for locating warehouses for multiple products. They also extended the model previously developed by Jayaraman for predetermined plant locations. Amiri (2006) formulated a model for determining the optimal numbers, locations, and capacities of singleproduct plants and warehouses. In a review study by Bektas (2006), it was stated that the features of the traveling salesman problem make it perfect for modeling many reallife issues and that it can be extended to a wide range of vehicle routing problems by adding a few auxiliary constraints. This study also provided an overview of the challenges in modeling a large number of traveling salesmen.
In a study by Koskosidis et al. (1992) they integrated vehicle routing and scheduling problems, stating that integrated models are better fitted to represent a set of problems. They presented a classification method for different types of these problems and examined the effectiveness of exact and heuristic methods in solving those (Taknezhad and Vahedi, 2021; Scott, 2010). Javid and Azad (2010) developed a locationinventoryrouting model for a supply chain with stochastic parameters, which does not use approximation. In a study by Shen et al. (2003), they showed that the approximation approach is only usable under certain assumptions and formulated their model as a convex mixedinteger program. Xu et al. (2013) developed a model for assigning retailers to distribution centers and planning their inventory levels in a supply chain where demand is stochastic and production centers are capacitated. In this model, inventory costs were also considered in allocation decisions.
Tancrez et al. (2012) proposed a model for a threelevel singlesource supply chain consisting of suppliers, distribution centers, and retailers that takes into account confidence and risk levels. In this model, waiting times are fixed and the objective is to determine the optimal number of suppliers and distribution centers, assign them to predetermined locations, and assign retailers to distribution centers. Sajjadi and Gharaei (2011) developed a model for multiproduct locationrouting in a supply chain with stochastic demand. In a study by Araghi et al. (2021), they presented a threelevel location model for a supply chain with Poisson demand. Hamedani et al. (2013) developed a multiobjective inventorylocation model that takes the risk level into account. Pahlevan et al. (2021) developed a location model for relocation planning in the aluminum industry and proposed several meta heuristic algorithms for solving this problem.
This paper examines the demand and capacity constraints that can be formulated when there is uncertainty at each level of a location problem. The paper also aims to provide a more extensive model of the locationinventory routing problem than those provided in previous works. In the proposed model, routing is treated as a dynamic problem taking into account the maintenance cost and the ordering costs at the end of each period. The model also includes budget constraints and a cap on the number of warehouses that can be established. To make the model more realistic, it is assumed that there is a degree of uncertainty in the demand parameters and the capacity of warehouses and retailers. In the end, the paper provides a robust optimization method for solving the problem under the defined uncertainty.
In the following, Section 2 defines the modeled problem, describes the modeling approach and the considered uncertainties, provides a method for solving the model, presents and discusses the results of applying the model on a case study, and finally presents the results of a sensitivity analysis. Section 3 present the conclusions and offers a few suggestions for future research.
2. MATHEMATICAL MODEL
The proposed model aims to find the best location of warehouses and the best routes between warehouses and stores in a supply chain consisting of a set of manufacturing plants, warehouses, and stores. The model has three levels as shown in Figure 1:
In the proposed model, it is assumed that there is no limit to the number of vehicles, the weight and volume capacity of each vehicle on each route is given, the location of stores is predetermined, the capacity of warehouses is limited, and each route starts from one warehouse and ends at the same warehouse and cannot go through another warehouse. The model also assumes that budget constraints cannot be violated and split delivery is not allowed (i.e. each demand is met by only one warehouse/vehicle). The proposed model is multiperiod and multiproduct and considered demand (at stores) to be stochastic. The objective of the model is to optimize the number of distribution centers and the routes between warehouses and stores for minimizing the cost of establishing warehouses, transportation costs, ordering and purchasing costs (for warehouses and stores), and maintenance costs.
2.1 Mathematical Modeling
The notations and symbols used in the model formulation are described below:
Sets and indices:

m: Index of manufacturing plants m = {1,2,…,M}

w: Index of warehouses w = {m+1, m+2, …, m+W}

s: Index of stores s = {m+w+1, m+w+2, …, m +w+S}

t: Index of periodst = {1, 2, …, T}

P: Index of products p = {1, 2, …, P}

i: Index of vehicles i = {1, 2, …, I}

k: Index of network nodes (warehouses and stores) $K=\left\{s{{\displaystyle \cup}}^{\text{}}w\right\}$
Parameters:

f_{w}: Fixed cost of establishing warehouse w

b_{wpt}: Maximum storage capacity of warehouse w for product p in period t

C_{wspt}: Cost of purchasing product p from manufacturing plant m by warehouse w in period t

d_{sw}: Distance between warehouse w and store s

A_{wpt}: Cost of ordering product p in period t for warehouse w

H_{pt}: Inventory cost of product p in period t for warehouse w

C_{it}: Cost of transport by vehicle i in period t (per unit distance traveled)

β_{i}: Capacity of vehicle i

D'_{sptv} Demand for product p at store s in period t (uncertain)

B: Maximum available budget

M*: A large number

ϑ: Discount coefficient

N: Number of warehouses to be established
Variables:

I_{wpt}: Inventory of product p at warehouse w at the end of period t

asw: =1 if store s is assigned to warehouse w, = 0 otherwise

X_{wspt}: The quantity of product p transferred to store s from warehouse w in period t (uncertain)

X'_{mwpt}: The quantity of product p transferred from manufacturing plant m to warehouse w in period t (uncertain)

Y_{w}: =1 if warehouse w is established, = 0 otherwise

V_{wspti}: =1 if vehicle i takes product p from warehouse w to store s in period t, =0 otherwise
Using the above notations, the proposed model is formulated as follows:
Equation (1) is the objective function, which consists of five cost terms: cost of purchasing product p by warehouse, cost of purchasing product p by stores, routing costs, ordering cost, maintenance cost, and fixed cost of establishing warehouses. Equation (2) is the storage capacity constraint. Equation (3) is the demand constraint. Equation (4) balances inventory levels in period t. Equation (5) limits the capacity of each warehouse. Equation (6) is the budget constraint. Equation (7) ensures that the number of stores does not exceed N because of budget constraints. Equation (8) is the vehicle capacity constrain. Equation (9) ensures that each store is assigned to only one warehouse. Equations (10), (11), (12) states that a store can have inventory at the end of the period and receive and send goods only if it is established. Equation (13) ensures that each store is assigned to one route only. Equation (14) states that if a vehicle is assigned to a warehouse, its first destination must be the first store assigned to that warehouse and its last destination should be the same warehouse. Equation (16) ensures that there is at least one route between each warehouse and its stores and from each store to the next. Equation (16) states that a store can be assigned to a warehouse only if an active route is passing through that store. Equation (17) ensures that each vehicle departs from its warehouse only once. Equation (18) defines the domains of nonnegative and binary variables.
2.2 Robust Optimization
In many mathematical programming (optimization) problems, there is a degree of uncertainty about the value of some parameters (e.g. demand, costs, operating conditions, etc.). In general, such uncertainties are addressed by the use of stochastic programming, fuzzy programming, robust programming, or combinations of these methods. For some problems, the provided solution must remain feasible for all or at least the majority of possible scenarios. In this study, we use the robust optimization approach of BenTal et al. (2009), which is an extension of the strict method of Sabela (2020). Therefore, we consider a bounded uncertainty interval for each certainty parameter and try to generate robust solutions that remain feasible for most possible cases of data uncertainty in the defined intervals. Assume that data is one of the uncertain parameters of the proposed model. For each data, we define an uncertainty interval as $data\in [dat{a}^{L}+dat{a}^{U}]$, based on which the nominal value $data=\frac{dat{a}^{L}+dat{a}^{U}}{2}$ and the deviation term $data=dat{a}^{U}=data$ of this uncertain parameter can be obtained. After determining the nominal values and deviation terms of all uncertain parameters, the robust counterpart of the model can be obtained.
2.3 Solution Method
2.3.1 Exact method
For the exact solution, an instance of the problem with one manufacturing plant, three potential warehouses, eight stores, one period, and one product was solved as a mixedinteger linear program with the GAMS software using the Baron method. The results of this method are provided in Tables 1 and 2. Among the potential warehouse locations, an optimal place was chosen as the place where the warehouse should be constructed. The demand estimated for each store, the quantity of goods received from the manufacturing plant, the routes assigned to the stores, and the total cost are also given in the Tables 1 and 2.
The best route found by the software is optimal13642875optimal, resulting in a total cost of 1.4360075 USD with a discount coefficient of 0.001.
2.3.2 Metaheuristic Method
Genetic Algorithm (GA) is a powerful search algorithm developed by Holland, which takes inspiration from the process of natural evolution. This evolutionary process begins with the generation of an initial population usually through random methods. The chromosomes within this population that meet certain conditions are subjected to crossover and mutation operators to produce new chromosomes. New chromosomes (offspring) are placed in the same pool as the original (parent) chromosomes. Since this pool is larger than the allowed population size, a number of chromosomes in the pool are picked by a selection operator to form the next generation. Usually, the new population is fitter than the previous one, or in other words, the population evolves from generation to generation. The search ends when either convergence or stopping conditions are met.
Control parameters of GA include population size (nPop), maximum iteration (MaxIt), crossover probability (pc), and mutation probability (pm). To determine the best values of these parameters for our problem, first, multiple levels of values that can be assigned to these parameters were chosen based on the choices made in previous studies. Then, Taguchi’s experimental design method was used to determine an appropriate experimental design (L9 (34)) for testing the effect of algorithm parameters on the solutions (Hartani et al., 2021;Koloba, 2020). Using this design, the effects were tested 27 times in MATLAB. Ultimately, the optimal parameter values were obtained in three levels as shown in Figure 2.
Since the goal was to minimize deviation from the objective function, it was decided to use the smallerthebetter variant of the Taguchi method. Using this method, it was determined that the maximum iteration should be set to the intermediate level, the population size should be set to the intermediate level, and the crossover and mutation probabilities should be set to their low level. The model was then solved with a GA coded in MATLAB with the parameter setting determined by the design of experiments. The model was also solved another time with the GAMS software. Table 3 presents the results obtained from MATLAB and GAMS, including the solution time in seconds and the total cost in million dollars.
2.4 Sensitivity Analysis
A sensitivity analysis was conducted to identify the parameters that strongly affect the model. This analysis was performed for the discount coefficient and fixed costs. For this purpose, the effect of a change in the discount coefficient on the objective function when other parameters are kept constant was investigated.
As shown in Figure 3, the choice of discount coefficient has a direct impact on the value of the objective function, and the total cost also decreases (increases) with the decrease (increase) in this parameter.
Figure 4 shows the simultaneous effect of giving different values to fixed costs and the discount coefficient on the objective function. It can be seen that the total cost decreases (increases) with the decrease (increase) in fixed costs and the discount coefficient. These results suggest that the designed mathematical model can find a good balance between fixed costs and variable costs.
3. CONCLUSIONS AND SUGGESTIONS FOR FUTURE RESEARCH
Treating all parameters of a supply chain model as deterministic variables will make it unrealistic and extremely vulnerable to any unforeseen change in parameter values. Therefore, to reach a more realistic locationrouting model, this study considered some parameters of the model to be inherently uncertain. Since locationrouting problems belong to the family of NPHard problems and become increasingly more complex with the rising number of nodes, which makes exact solution methods increasingly ineffective, the paper also provided a metaheuristic method for solving the model. The results showed little difference between the solutions of exact and metaheuristic methods for medium and smallsized instances of the problem.
In future studies, it is recommended to define different service levels for retailers based on their priority for the supply chain. Another interesting option is to expand the supply chain network into more levels and introduce suppliers and satellite or intermediate warehouses to the formulation.