# Providing a Mathematical Routing-Inventory Model for the Drug Supply Chain Considering the Travel Time Dependence and Perishability on Multiple Graphs

Authors: Madad Ali; Martha Ng; Rui Dias; Iskandar Muda; Riyadh Al-Obaidi; Barno Abdullaeva; Himanshu Sharma; Hussein Mohammed Esmail Abu Al-Rejal; Ali Thaeer Hammid
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

High turnover in the pharmaceutical industry, the location of placebos in urban spaces, and the high rate of corruption of products in this industry are the distinguishing features of the drug supply chain. Thus, to survive and maintain competitive advantages in the current business environment, managers are active in this background to implement the theoretical foundations of supply chain management. One of the influential areas in this category is integrated inventory management, inventory control, and vehicle routing. Therefore, this study mainly aims to analyze and define a routing problem, inventory in the drug supply chain with perishable products, and travel time dependence on multiple graphs with travel time dependence. The Box–Jenkins forecasting method has been utilized to meet the study’s aim and deal with demand uncertainty. This method can identify the best pattern governing the data. Finally, the mathematical model is validated, and managerial perspectives are provided. The study results demonstrated the possibility of achieving cost-saving and reducing product spoilage. Applying the solutions of this model can provide some inherent social and environmental advantages, including reducing traffic load and emissions.


---

## Body

## 1. Introduction

In a world of growing healthcare systems with massive economies, it is essential to have some means of controlling and reducing expenses. Global healthcare expenditures are expected to rise by approximately 5.4%, from 7.724 trillion dollars in 2017 to 10.059 trillion dollars in 2022 [1, 2]. Currently, pharmaceutical products account for 20–30% of the total cost of healthcare systems [3]. This highlights the amount of cost-saving that can be achieved by reducing pharmaceutical product expenditures. According to Aptel and Pourjalili [4], managing pharmaceutical supply chains is one of the critical managerial tasks in the pharmaceutical industry. It has been suggested that integrating supply chain decisions can help achieve significant cost savings in this area.In addition to their substantial economic size, pharmaceutical supply chains are characterized by the complexity of supply chain planning and control. Factors contributing to this complexity include production and distribution requirements such as good manufacturing practice (GMP) and good distribution practice (GDP), a considerable number of customers, and the dispersed distribution of retailers. Good manufacturing practice (GMP) is regarded as a system to ensure that products are consistently controlled and produced based on quality measures and criteria. It is mainly developed to reduce the risks involved in pharmaceutical production that cannot be easily eradicated via experimenting with the final product [5].According to Kapoor et al. [5], other challenges of pharmaceutical supply chains include the difficulty of coordination between chain components, shortages, the expiration of products, and the requirements for safe distribution. These characteristics distinguish pharmaceutical supply chains from many other chains. Efficient and effective management of pharmaceutical supply chains can boost the development of almost all business processes in these chains. Therefore, it is crucial to use the principles of supply chain management in this industry for cost-saving and profitability improvement purposes. Essentially, a pharmaceutical supply chain is the integration of all activities related to the flow and transformation of raw materials into pharmaceutical products and their delivery to end customers, as well as the sharing and flow of information among chain components to achieve sustainable competitive advantage [6–8]. In general, pharmaceutical supply chain management problems can be divided into two domains of resource management and process development or improvement [9–11]. In a study by Mousazadeh et al. [8], they designed a four-level supply chain model for making strategic decisions regarding the construction of production centers and tactical decisions regarding the flow of materials in the medium-term planning horizon. Shang et al. [9] provided a solution for redesigning a pharmaceutical company’s supply chain network to reduce costs and increase customer satisfaction. Marques et al. [10] modeled the design and manufacturing process of a new pharmaceutical product with an integrated linear programming approach. In a study by Zahiri et al. [11], they designed a resilient pharmaceutical supply chain network. Most articles in this field deal with inventory management problems at the tactical level. For example, Candan and Yazgan [12] formulated the pharmaceutical supply chain inventory problem as an integrated linear program to maximize net profit. These researchers used classical and vendor-managed inventory (VMI) problem instances to demonstrate the effectiveness of their model. Uthayakumar and Priyan [13] developed a model for integrated production-inventory planning in a two-level supply chain comprised of factories and hospitals. Jillson et al. [14] formulated a model for dynamic routing of high-value, low-volume products at the operational level. In general, articles on decision-making at the tactical level (e.g., inventory management) and at the operational level (e.g., routing) separately or together constitute a small portion of the subject literature [15, 16].Many assumptions about classical routing problems no longer apply to the world’s major cities and metropolitan areas. For example, many modern urban transportation networks are designed to have multiple viable routes between any two nodes. Also, since the traffic load depends on the time of day, the travel time will likely depend on when the trip occurs. Thus, unlike in classic routing problems, it is possible to travel longer distances in a shorter time at certain hours of the day. The problem of vehicle routing with time-dependent travels was first introduced in 1992 by Malandraki and Daskin [15], who formulated the traveling salesman problem, which is a special case of the routing problem, as a complex integer model. In their proposed model, a time window can also be considered. The main drawback of their work was that travel time was considered a stepwise function of the time of day in which the travel takes place, which can lead to the violation of the first-in-first-out (FIFO) principle. To resolve this issue, they suggested considering the notion of vehicles waiting at the nodes, which can be costly and impractical. In 2003, Ichoua et al. [16] proposed a conversion-based method to address this issue. Their method is an algorithm that converts the travel speed function into a time-dependent travel time function. In a study by Setak et al. [17], they formulated the concept of multiple paths between nodes in urban transportation networks as a multi-graph. Later, Huang et al. [18] and Patoghi et al. [19], and Tikani and Setak [20] considered this issue in their articles. Ignoring the time dependence of travel time, especially when this dependence is strong, can significantly undermine the validity, applicability, and optimality of the obtained solutions. While recent years have seen increasing interest in this domain of ​​research in order to make vehicle routing models more realistic, there is still much work to be done in this area.Poor inventory management in pharmaceutical supply chains can increase the spoilage rate and, therefore, the supply chain costs. In North America, for example, drug spoilage cost the pharmaceutical industry about 396 million dollars in 2018, and this cost is projected to increase at a rate of 5.4% until 2025. Considering the decreasing profit margin of pharmaceutical supply chains, it is imperative to consider these products as perishable items during planning to reduce inventory waste costs. Only a small number of articles in the pharmaceutical supply chain management field have considered inventories perishable. These include the works of Papageorgiou et al. [21] on product development and capacity planning and Shen et al. [22] on inventory management. In the field of vehicle routing, the studies that have considered perishable inventories include the studies of Le et al. [23], Coelho and Laporte [24], and Rahimi et al. [25]. Therefore, it is important to consider products as perishable items in future formulations.One of the main assumptions of nonstochastic (deterministic) vehicle routing and inventory problems is that all parameters and data are given or obtained at certain values [26]. This assumption has a significant drawback in that the solutions become infeasible or suboptimal when the real-world values of parameters fall outside their assumed ranges. In the literature on inventory routing for perishable products, three methods have been used to deal with uncertainty: stochastic modeling, robust modeling, and fuzzy modeling. For example, Goli and Malmir [27], and Goli and Mohammadi [28], used stochastic methods to control uncertainties in their plans. Goli et al. [29] used a robust method to control the uncertainty in demand. In an article by Rahimi et al. [25], the uncertainty was controlled by a fuzzy approach, specifically by using a probabilistic fuzzy triangular function with three sets of optimistic, pessimistic, and most probable values. As can be seen, a growing number of researchers working in the field have shown an interest in modeling uncertainty in their models to make their outputs more realistic. Nevertheless, given the varying circumstances of different problems, there is still much work to be done in this area [30, 31]. Also, the studies of Shaabani and Kamalabadi [32], Sepehri and Sazvar [33], Soysal et al. [34], Jafarkhan and Yaghoubi [35], and Crama et al. [36] have proposed new methods for dealing with uncertainty based on fuzzy numbers.Reviewing the related studies and sources, it can be inferred that inventory routing problems are inherently np-hard, and uncertainties should be controlled without losing attention to the complexity of the model. Hence, considering the abovementioned issues surrounding pharmaceutical supply chains, to fill the existing gaps in the literature, this study attempts to provide an integrated inventory-routing model for perishable pharmaceutical products with uncertain demand and time-dependent demand travel times. The main contributions of the study to the respective field can be summarized as follows:(1)
The study utilizes integrated pharmaceutical supply chain decisions, which can substantially reduce the cost of the entire chain.(2)
The time dependence of travel time is considered on a multi-graph with a heterogeneous transport fleet, which can assist in making the model more realistic and more applicable to real-world problems.(3)
A spoilage rate for the inventory is determined, which can help reduce inventory waste costs.(4)
The Box-Jenkins forecasting methodology is taken into account to predict the uncertain demand in a case study. This method estimates the model that best explains the data while respecting the principle of parsimony (Occam’s razor). Furthermore, given that the formulations are applied to a case study, the available data are employed to prevent the problem solving, which is NP-Hard, from becoming more complex and time-consuming.


## 2. Materials and Methods

In this section, the main assumptions of the model presented in the third section are described, which take inspiration from the real-world situation of the considered case. Basically, the proposed model is based on the pharmaceutical distribution companies’ inventory and distribution management processes. This chain has two distinctive features: the perishability of products and the dispersed distribution of customers over the entire urban transportation network. The second feature, which is supposed to facilitate the accessibility of customers, i.e., pharmacies, clinics, and hospitals, leads to the dependence of travel time on the time of day when the travel takes place and also the superiority of a multiroute approach. Of course, the mentioned features do not cause the model to lose generality, and it can still be adapted to other chains such as cold supply chains and organ transplant supply chains.Generally, the warehouse must have sufficient inventory to respond to the requested orders. This inventory level should be maintained according to projected demand and the spoilage rate of products. Orders are loaded onto the heterogeneous transport fleet at the beginning of the delivery horizon. The route sequences and active routes for loaded vehicles are determined according to the traffic load around customers during the day when the trip takes place. The vehicles return to the warehouse after serving the last customer.This problem is formulated on a complete multigraph G(N, A), where N represents the set of nodes and A represents the set of edges connecting the nodes. Each edge is denoted by (i, j, p), representing the p-th path between the nodes i and j. Node 0 represents the supplier, and the remaining nodes represent the customers. The set of vehicles is denoted by vt∈1,2,…,V. Each type of vehicle has a certain capacity, denoted by capvt. The number of vehicles in the warehouse is denoted by kvt. Each vehicle can take only a single route in each time period of the planning horizon. Also, each customer can be served by only one vehicle in each time period of the planning horizon. The service time is denoted by si. All activated vehicles start their trips from the warehouse simultaneously at the beginning of the delivery horizon and return to the warehouse at the end of their service. An inventory holding cost of hi is incurred for each unit of product that is stored in each node in each period. Since each product has a limited shelf life, wri percent of products spoil in each period. For each unit of spoiled product, the system incurs a cost denoted by pr, which can be as high as the product's unit price.Travel time on each edge depends on the time of day when the trip takes place. So as to satisfy the FIFO principle, the travel speed function was converted to travel time. The main reason for this action is that the mentioned conversion makes the model more realistic and accurate, given the fact that the change in travel speed is more tangible and substantial.The purpose of the model is to determine the optimal delivery routes, the optimal quantity of goods to be delivered to each customer, and the optimal inventory level of each node in each planning horizon period such that routing, inventory, and spoilage costs are minimized. It is assumed that routing costs are comprised of a variable time-dependent cost and a fixed cost for using the vehicle.


### 2.1. Assumptions

To carry out this study, it is assumed that the entire demand in each period must be fully met. Having historical demand data, the Box–Jenkins forecasting method is utilized to predict the demand in the future periods of the planning horizon for use in the model table1.Table 1
Indices, variables, and parameters of the proposed model.
Indices and setsk,j,iIndex for nodes; node 0 represents the central warehouse and other nodes represent customers.PtIndex of periods of the planning horizonvtIndex of vehicle typesPIndex of the p-th path between two nodesrThe number of intervals into which the delivery time horizon is divided such that the travel speed on the route remains constant in each interval.mIndex of intervals in which travel time is a piecewise linear function with a fixed slope in each interval.Parametersbi,j,vtp,m,bi,j,vtp,m+1The m-th interval of the travel time functionti,jp,r,ti,jp,r+1The r-th interval of the travel speed functionvi,j,vtp,rTravel speed on pathp between nodes i and j for vehicle vt in interval rdi,jpThe distance between nodes i andj via path pdepi,j,vtp,rThe time that vehicle vt needs to depart from node i toward nodej in order to arrive at its destination at a specified time via path p in interval rarivi,j,vtp,mThe time that vehicle vt arrives at nodej via path p if it departs from node i in interval m.trti,j,vtp,mTravel time from node i to nodej via path p for vehicle vt if it departs from node i in interval mθi,j,vtp,mThe slope of the piecewise linear function of travel time for the node pair (i, j), vehicle vt, pathp, and interval mηi,j,vtp,mThe intercept of the piecewise linear function of travel time for the node pair (i, j), vehicle vt, pathp, and interval mSiService time for customer idemiptThe demand of customer i in period ptcapvtThe capacity of a vehicle vtcapiWarehousing capacity of node ifcvtFixed cost of using vehicle vtVCvtThe time-dependent variable of using vehicle vthiCost of storing the product in node i for one periodPrThe unit price of the productwriProduct spoilage rate when stored in node iiiliInitial inventory level of node i at the beginning of the planning horizon (it is assumed that the initial inventory level of the central warehouse is 0).kvtNumber of vehicles of type vtL,ULength of the delivery time horizon in each period of the planning horizonVariablesZObjective function valuexi,jvt,pt=1 if the edge connecting node i to nodej is selected as the optimal route for vehicle vt in period pt; = 0 otherwisexi,j,vtp,pt=1 if the edge connecting node i to nodej via path p is selected as the optimal route for vehicle vt in period pt;=0 otherwisexi,j,vtp,m,pt=1 if the edge connecting node i to nodej via path p is selected as the optimal route for vehicle vt in interval m in period pt; = 0 otherwisedti,j,vtp,m,ptTime of departure of vehicle vt from node i toward nodej via path p in interval m in period ptfi,jptThe amount of product transported on the activated edge connecting node i to nodej in period ptfi,jp,ptThe amount of product transported on the activated edge connecting node i to nodej via path p in period ptfi,j,vtp,ptThe amount of product transported by vehicle vt on the activated edge connecting node i to nodej via path p in period ptiliptInventory level of node i at the end of period ptwasiptThe amount of product spoiled in the warehouse of node i during the period ptqiptThe amount of product delivered to customer i in period ptrpdptThe amount of product arriving at the central warehouse during the period pt


### 2.2. Objective Function

The objective function of the model is to minimize costs. These costs are divided into two categories: routing costs and inventory costs. Routing costs include the fixed cost of using the vehicles, which depends on the type of vehicle, and a variable cost, which depends on the number of hours the vehicle is to be used over the course of the planning horizon. Inventory costs consist of inventory holding costs and product spoilage costs. Therefore, the objective function is formulated as follows:(1)minZ=∑vt=1Vfcvt.∑i=0N∑j=0i≠jN∑m=1M∑pt=1T∑p=1Pxi,j,vtp,m,pt+∑vt=1VVCvt.∑i=1N∑p=1P∑m=1M∑pt=1Tθi,0,vtp,mdti,0,vtp,m,pt+ηi,j,vtp,m.xi,0,vtp,m,pt+∑i=0Nhi.∑pt=1Tilipt+pr.∑i=0N∑pt=1Twasipt.


### 2.3. Model Constraints

The objective function formulated in the previous subsection is based on the following constraints:(2)il0pt=il0pt−1+rpdpt−was0pt−∑j=1Nqjpt∀pt∈T1,(3)il01=rpd1−was01−∑i=1Nqi1,(4)il0pt=il0pt−1+qipt−was0pt−demipt∀pt∈T1,i∈N0,(5)ili1=iili+qipt−wasipt−demipt∀i∈N0,(6)ilipt≤capi∀i∈N,pt∈T,(7)wasipt=wri.ilipt∀i∈N,pt∈T,(8)rpdpt=∑i=0Nwasipt+∑i=1Ndemipt∀pt∈T,(9)∑j=1Nx0,jvt,pt≤kvt∀vt∈V,pt∈T,(10)∑i=0N∑vt=1Vxi,jvt,pt=1∀j∈N0,pt∈T,i≠j,(11)∑j=0N∑vt=1Vxi,jvt,pt=1∀i∈N0,pt∈T,i≠j,(12)∑i=0Nxi,jvt,pt=∑k=0Nxj,kvt,pt∀j∈N0,vt∈V,pt≠T,(13)∑p=1Pxi,j,vtp,pt=xi,jvt,pt∀i,j∈N,vt∈V,pt∈T,(14)∑m=1Mxi,j,vtp,m,pt=xi,j,vtp,pt∀i,j∈N,i≠j,p∈P,vt∈V,pt∈T,(15)xi,i,vtp,m,pt=0∀i∈N,vt∈V,p∈P,m∈M,pt∈T,(16)∑i=0Nfi,jpt−∑k=0Nfj,kpt=qjpt∀j∈N0,pt∈T,(17)∑p=1P∑vt=1Vfi,j,vtp,pt=fi,jpt∀i,j∈N,i≠j,pt∈T,(18)qjpt≤fi,j,vvp,pt≤capvt−qipt∀i,j∈N,i≠j,p∈P,vt∈V,pt∈T(19)dt0,j,vtp,m,pt=L.x0,j,vtp,m,pt∀j∈N0,vt∈V,p∈P,m∈M,pt∈T,(20)bi,j,vtp,m.xi,j,vtp,m,pt≤dti,j,vtp,m,pt≤bi,j,vtp,m+1.xi,j,vtp,m,pt∀i,j∈N,i≠j,vt∈V,p∈P,m∈M,pt∈T,(21)∑i=0N∑p=1P∑m=1M∑vt=1Vθi,j,vtp,m.dti,j,vtp,m,pt+ηi,j,vtp,m.xi,j,vtp,m,pt+∑i=0N∑p=1P∑m=1M∑vt=1Vdti,j,vtp,m,pt+∑i=0N∑p=1P∑m=1M∑vt=1Vdtj,k,vtp,m,pt∀j∈N0,j≠i,j≠k,pt∈T,(22)∑p=1P∑m=1M∑vt=1Vθi,o,vtp,m.dti,o,vtp,m,pt+ηi,o,vtp,m.xi,o,vtp,m,pt≤U∀i∈N0pt∈T,(23)Z,dti,j,vtp,m,pt,fi,jpt,fi,jp,pt,fi,j,vtp,pt,ilipt,wasipt,qipt,rpdpt≥0,(24)xi,jvt,pt,xi,j,vtp,pt,xi,j,vtp,m,pt=oorl.Constraint (1) maintains the inventory balance of the central warehouse. Constraint (2) computes the inventory level at the end of the first period based on the assumption that the inventory level at the beginning of the planning horizon is zero. Constraint (3) maintains the inventory balance of customer warehouses. Constraint (4) determines the inventory level of customers at the end of the first period, assuming that their inventory is empty at the beginning of the planning horizon. Constraint (5) limits the capacity of warehouses. Constraint (6) calculates the amount of product spoiled in warehouses. Constraint (7) determines the amount of product received in the central warehouse in each period. Constraint (8) limits the number of activated vehicles of any type to the number of vehicles of that type that are available. Constraints (9) and (10) ensure that each node is visited exactly once. Constraint (11) forces activated vehicles to depart toward another customer or the central warehouse upon serving a customer. Constraint (12) states that if a vehicle of type vt is selected to travel between two nodes (i, j), it will travel only one route (path p) between the two nodes. Constraint (13) states that if a path p between two nodes (i, j) is activated for a vehicle of type vt, this travel must take place exactly in the m-th time interval bi,j,vtp,m,bi,j,vtp,m+1. Constraint (14) prevents the formation of loops. Constraint (15) balances the flow of products between nodes to satisfy customer demand. Constraint (16) ensures that products are transported on a single route between (i, j). Constraint (17) makes sure that the capacity of each type of vehicle is respected while serving customers. Constraint (18) states that no vehicle can leave the central warehouse before the start of the delivery time horizon. Constraint (19) limits the time of departure from node i toward node j on path p to the m-th time interval specified in Constraint (20). In Constraint (21), the time of departure from each node is computed by summing the time of departure from the previous node with the travel time between the two nodes and the service time of that node. Constraint (22) states that the delivery to customers must be completed before the end of the delivery time horizon. Constraints (23) and (24) specify the type of variables.


### 2.4. Solution Procedure

The proposed routing-inventory model has complexities in various respects. On the one hand, the dependence of travel time on the time of day and the path between the nodes is available. This issue can be resolved by utilizing a stepwise time-dependent travel time function. Nonetheless, applying this approach results in a violation of the FIFO principle. On the other hand, we have to cope with uncertainty in the demand parameter without adding to the complexity of the model. To resolve these issues, a two-step approach was proposed. In the first step, the time-dependent travel time function was computed while observing the FIFO principle for each route. The output of this step will be a time-dependent routing-inventory model. In the next step, the Box–Jenkins forecasting method was employed to predict the demand. Adopting this approach prevents the model from getting more complex while observing the principle of parsimony. The outputs of this procedure will be the inputs of the model provided in Section3.


### 2.5. Time-dependent Travel Time Function Calculations

Considering a constant length for time-dependent travel time leads to the violation of the FIFO principle, which states that if vehicle A departs from node i toward nodej earlier than an identical vehicle like B does the same, vehicle A must reach its destination sooner. In an article by Ichoua et al. [16], they presented a technique for solving this problem based on the decomposition of travel speed over consecutive periods for each day. In this way, they divided the day into multiple time periods over which travel speed remains constant. This approach gives a time-dependent travel speed function in the form of a stepwise function. Thus, travel time can be calculated from the travel speed function with an algorithm based on the decomposition of travel speed over consecutive periods as described below.It is assumed that the length of the delivery time horizon [L, U] is divided intor intervals with endpoints T=t1,t2,…,tr, where t1=L and tr+1=u, such that speed over the interval r remains constant at vr. Since there can be more than one route between every two nodes, the travel time function corresponding to the travel speed function for each route will be a piecewise linear function.


### 2.6. Travel Speed Conversion

The main idea of this method is that vehicles do not travel the entire route between two nodes at a constant speed but instead change their speed at two or more time intervals. For example, a vehicle that starts its trip from node i to nodej at the time t0 of the day may travel one part of the route d1 at speed v1 over time t1−t0, travel another part of the route d2 at speed v2 over time t2−t1, and travel the remainder of the route d3=d−d2+d1 at speed v3.


### 2.7. Travel Time Calculation with the Travel Speed Function

For a vehicle leaving node i at a timet0∈Tr=ti,jp,r,ti,jp,r+1 to reach node j via path p, the travel time is calculated through the following procedure:1-1- Sett to t0.1-2- Setd to dij.1-3-Sett′ to t+d/vr.2-Perform the following loop as long ast′>tijpr.2-1-Setd to d−vrtijp,r+1−t.2-2-Sett to tijp,r+1.2-3-Sett′ to t+d/vr+1.2-4-Setr to r+1.3-Travel time =t′−t0.


### 2.8. Box-Jenkins Forecasting Method

A time series is an ordered sequence of data over time. One inherent characteristic of time series is the dependence or correlation of observations, which makes the order of observations quite important. Hence, statistical methods that operate based on the assumption that data are independent over time are not applicable to time series. The statistical methodology developed for working with and interpreting time series is called time series analysis. Typically, the goal of time series analysis is to describe, model, or predict data or examine the dynamics of a variable over time. In the current study, the uncertainty in demand is modeled utilizing the Box–Jenkins forecasting methodology, involving fitting an auto-regressive integrated moving average (ARIMA) model [37–39]. An ARIMA model itself is characterized by its autoregressive and moving average processes and a degree of integration. The moving average process models the observed time series Zt as a linear combination of a sequence of uncorrelated random variables. The autoregressive process models the present value of the time series as a function of its previous values plus a random effect. Integration is used for time series that do not have a fixed mean to make analyses more reliable. Therefore, if a time series is stationary after d times of first-order differencing and is modeled as ARMA (p, q), then the original time series is an autoregressive integrated moving average time series in the form of ARIMA (p, d, q).


## 2.1. Assumptions

To carry out this study, it is assumed that the entire demand in each period must be fully met. Having historical demand data, the Box–Jenkins forecasting method is utilized to predict the demand in the future periods of the planning horizon for use in the model table1.Table 1
Indices, variables, and parameters of the proposed model.
Indices and setsk,j,iIndex for nodes; node 0 represents the central warehouse and other nodes represent customers.PtIndex of periods of the planning horizonvtIndex of vehicle typesPIndex of the p-th path between two nodesrThe number of intervals into which the delivery time horizon is divided such that the travel speed on the route remains constant in each interval.mIndex of intervals in which travel time is a piecewise linear function with a fixed slope in each interval.Parametersbi,j,vtp,m,bi,j,vtp,m+1The m-th interval of the travel time functionti,jp,r,ti,jp,r+1The r-th interval of the travel speed functionvi,j,vtp,rTravel speed on pathp between nodes i and j for vehicle vt in interval rdi,jpThe distance between nodes i andj via path pdepi,j,vtp,rThe time that vehicle vt needs to depart from node i toward nodej in order to arrive at its destination at a specified time via path p in interval rarivi,j,vtp,mThe time that vehicle vt arrives at nodej via path p if it departs from node i in interval m.trti,j,vtp,mTravel time from node i to nodej via path p for vehicle vt if it departs from node i in interval mθi,j,vtp,mThe slope of the piecewise linear function of travel time for the node pair (i, j), vehicle vt, pathp, and interval mηi,j,vtp,mThe intercept of the piecewise linear function of travel time for the node pair (i, j), vehicle vt, pathp, and interval mSiService time for customer idemiptThe demand of customer i in period ptcapvtThe capacity of a vehicle vtcapiWarehousing capacity of node ifcvtFixed cost of using vehicle vtVCvtThe time-dependent variable of using vehicle vthiCost of storing the product in node i for one periodPrThe unit price of the productwriProduct spoilage rate when stored in node iiiliInitial inventory level of node i at the beginning of the planning horizon (it is assumed that the initial inventory level of the central warehouse is 0).kvtNumber of vehicles of type vtL,ULength of the delivery time horizon in each period of the planning horizonVariablesZObjective function valuexi,jvt,pt=1 if the edge connecting node i to nodej is selected as the optimal route for vehicle vt in period pt; = 0 otherwisexi,j,vtp,pt=1 if the edge connecting node i to nodej via path p is selected as the optimal route for vehicle vt in period pt;=0 otherwisexi,j,vtp,m,pt=1 if the edge connecting node i to nodej via path p is selected as the optimal route for vehicle vt in interval m in period pt; = 0 otherwisedti,j,vtp,m,ptTime of departure of vehicle vt from node i toward nodej via path p in interval m in period ptfi,jptThe amount of product transported on the activated edge connecting node i to nodej in period ptfi,jp,ptThe amount of product transported on the activated edge connecting node i to nodej via path p in period ptfi,j,vtp,ptThe amount of product transported by vehicle vt on the activated edge connecting node i to nodej via path p in period ptiliptInventory level of node i at the end of period ptwasiptThe amount of product spoiled in the warehouse of node i during the period ptqiptThe amount of product delivered to customer i in period ptrpdptThe amount of product arriving at the central warehouse during the period pt


## 2.2. Objective Function

The objective function of the model is to minimize costs. These costs are divided into two categories: routing costs and inventory costs. Routing costs include the fixed cost of using the vehicles, which depends on the type of vehicle, and a variable cost, which depends on the number of hours the vehicle is to be used over the course of the planning horizon. Inventory costs consist of inventory holding costs and product spoilage costs. Therefore, the objective function is formulated as follows:(1)minZ=∑vt=1Vfcvt.∑i=0N∑j=0i≠jN∑m=1M∑pt=1T∑p=1Pxi,j,vtp,m,pt+∑vt=1VVCvt.∑i=1N∑p=1P∑m=1M∑pt=1Tθi,0,vtp,mdti,0,vtp,m,pt+ηi,j,vtp,m.xi,0,vtp,m,pt+∑i=0Nhi.∑pt=1Tilipt+pr.∑i=0N∑pt=1Twasipt.


## 2.3. Model Constraints

The objective function formulated in the previous subsection is based on the following constraints:(2)il0pt=il0pt−1+rpdpt−was0pt−∑j=1Nqjpt∀pt∈T1,(3)il01=rpd1−was01−∑i=1Nqi1,(4)il0pt=il0pt−1+qipt−was0pt−demipt∀pt∈T1,i∈N0,(5)ili1=iili+qipt−wasipt−demipt∀i∈N0,(6)ilipt≤capi∀i∈N,pt∈T,(7)wasipt=wri.ilipt∀i∈N,pt∈T,(8)rpdpt=∑i=0Nwasipt+∑i=1Ndemipt∀pt∈T,(9)∑j=1Nx0,jvt,pt≤kvt∀vt∈V,pt∈T,(10)∑i=0N∑vt=1Vxi,jvt,pt=1∀j∈N0,pt∈T,i≠j,(11)∑j=0N∑vt=1Vxi,jvt,pt=1∀i∈N0,pt∈T,i≠j,(12)∑i=0Nxi,jvt,pt=∑k=0Nxj,kvt,pt∀j∈N0,vt∈V,pt≠T,(13)∑p=1Pxi,j,vtp,pt=xi,jvt,pt∀i,j∈N,vt∈V,pt∈T,(14)∑m=1Mxi,j,vtp,m,pt=xi,j,vtp,pt∀i,j∈N,i≠j,p∈P,vt∈V,pt∈T,(15)xi,i,vtp,m,pt=0∀i∈N,vt∈V,p∈P,m∈M,pt∈T,(16)∑i=0Nfi,jpt−∑k=0Nfj,kpt=qjpt∀j∈N0,pt∈T,(17)∑p=1P∑vt=1Vfi,j,vtp,pt=fi,jpt∀i,j∈N,i≠j,pt∈T,(18)qjpt≤fi,j,vvp,pt≤capvt−qipt∀i,j∈N,i≠j,p∈P,vt∈V,pt∈T(19)dt0,j,vtp,m,pt=L.x0,j,vtp,m,pt∀j∈N0,vt∈V,p∈P,m∈M,pt∈T,(20)bi,j,vtp,m.xi,j,vtp,m,pt≤dti,j,vtp,m,pt≤bi,j,vtp,m+1.xi,j,vtp,m,pt∀i,j∈N,i≠j,vt∈V,p∈P,m∈M,pt∈T,(21)∑i=0N∑p=1P∑m=1M∑vt=1Vθi,j,vtp,m.dti,j,vtp,m,pt+ηi,j,vtp,m.xi,j,vtp,m,pt+∑i=0N∑p=1P∑m=1M∑vt=1Vdti,j,vtp,m,pt+∑i=0N∑p=1P∑m=1M∑vt=1Vdtj,k,vtp,m,pt∀j∈N0,j≠i,j≠k,pt∈T,(22)∑p=1P∑m=1M∑vt=1Vθi,o,vtp,m.dti,o,vtp,m,pt+ηi,o,vtp,m.xi,o,vtp,m,pt≤U∀i∈N0pt∈T,(23)Z,dti,j,vtp,m,pt,fi,jpt,fi,jp,pt,fi,j,vtp,pt,ilipt,wasipt,qipt,rpdpt≥0,(24)xi,jvt,pt,xi,j,vtp,pt,xi,j,vtp,m,pt=oorl.Constraint (1) maintains the inventory balance of the central warehouse. Constraint (2) computes the inventory level at the end of the first period based on the assumption that the inventory level at the beginning of the planning horizon is zero. Constraint (3) maintains the inventory balance of customer warehouses. Constraint (4) determines the inventory level of customers at the end of the first period, assuming that their inventory is empty at the beginning of the planning horizon. Constraint (5) limits the capacity of warehouses. Constraint (6) calculates the amount of product spoiled in warehouses. Constraint (7) determines the amount of product received in the central warehouse in each period. Constraint (8) limits the number of activated vehicles of any type to the number of vehicles of that type that are available. Constraints (9) and (10) ensure that each node is visited exactly once. Constraint (11) forces activated vehicles to depart toward another customer or the central warehouse upon serving a customer. Constraint (12) states that if a vehicle of type vt is selected to travel between two nodes (i, j), it will travel only one route (path p) between the two nodes. Constraint (13) states that if a path p between two nodes (i, j) is activated for a vehicle of type vt, this travel must take place exactly in the m-th time interval bi,j,vtp,m,bi,j,vtp,m+1. Constraint (14) prevents the formation of loops. Constraint (15) balances the flow of products between nodes to satisfy customer demand. Constraint (16) ensures that products are transported on a single route between (i, j). Constraint (17) makes sure that the capacity of each type of vehicle is respected while serving customers. Constraint (18) states that no vehicle can leave the central warehouse before the start of the delivery time horizon. Constraint (19) limits the time of departure from node i toward node j on path p to the m-th time interval specified in Constraint (20). In Constraint (21), the time of departure from each node is computed by summing the time of departure from the previous node with the travel time between the two nodes and the service time of that node. Constraint (22) states that the delivery to customers must be completed before the end of the delivery time horizon. Constraints (23) and (24) specify the type of variables.


## 2.4. Solution Procedure

The proposed routing-inventory model has complexities in various respects. On the one hand, the dependence of travel time on the time of day and the path between the nodes is available. This issue can be resolved by utilizing a stepwise time-dependent travel time function. Nonetheless, applying this approach results in a violation of the FIFO principle. On the other hand, we have to cope with uncertainty in the demand parameter without adding to the complexity of the model. To resolve these issues, a two-step approach was proposed. In the first step, the time-dependent travel time function was computed while observing the FIFO principle for each route. The output of this step will be a time-dependent routing-inventory model. In the next step, the Box–Jenkins forecasting method was employed to predict the demand. Adopting this approach prevents the model from getting more complex while observing the principle of parsimony. The outputs of this procedure will be the inputs of the model provided in Section3.


## 2.5. Time-dependent Travel Time Function Calculations

Considering a constant length for time-dependent travel time leads to the violation of the FIFO principle, which states that if vehicle A departs from node i toward nodej earlier than an identical vehicle like B does the same, vehicle A must reach its destination sooner. In an article by Ichoua et al. [16], they presented a technique for solving this problem based on the decomposition of travel speed over consecutive periods for each day. In this way, they divided the day into multiple time periods over which travel speed remains constant. This approach gives a time-dependent travel speed function in the form of a stepwise function. Thus, travel time can be calculated from the travel speed function with an algorithm based on the decomposition of travel speed over consecutive periods as described below.It is assumed that the length of the delivery time horizon [L, U] is divided intor intervals with endpoints T=t1,t2,…,tr, where t1=L and tr+1=u, such that speed over the interval r remains constant at vr. Since there can be more than one route between every two nodes, the travel time function corresponding to the travel speed function for each route will be a piecewise linear function.


## 2.6. Travel Speed Conversion

The main idea of this method is that vehicles do not travel the entire route between two nodes at a constant speed but instead change their speed at two or more time intervals. For example, a vehicle that starts its trip from node i to nodej at the time t0 of the day may travel one part of the route d1 at speed v1 over time t1−t0, travel another part of the route d2 at speed v2 over time t2−t1, and travel the remainder of the route d3=d−d2+d1 at speed v3.


## 2.7. Travel Time Calculation with the Travel Speed Function

For a vehicle leaving node i at a timet0∈Tr=ti,jp,r,ti,jp,r+1 to reach node j via path p, the travel time is calculated through the following procedure:1-1- Sett to t0.1-2- Setd to dij.1-3-Sett′ to t+d/vr.2-Perform the following loop as long ast′>tijpr.2-1-Setd to d−vrtijp,r+1−t.2-2-Sett to tijp,r+1.2-3-Sett′ to t+d/vr+1.2-4-Setr to r+1.3-Travel time =t′−t0.


## 2.8. Box-Jenkins Forecasting Method

A time series is an ordered sequence of data over time. One inherent characteristic of time series is the dependence or correlation of observations, which makes the order of observations quite important. Hence, statistical methods that operate based on the assumption that data are independent over time are not applicable to time series. The statistical methodology developed for working with and interpreting time series is called time series analysis. Typically, the goal of time series analysis is to describe, model, or predict data or examine the dynamics of a variable over time. In the current study, the uncertainty in demand is modeled utilizing the Box–Jenkins forecasting methodology, involving fitting an auto-regressive integrated moving average (ARIMA) model [37–39]. An ARIMA model itself is characterized by its autoregressive and moving average processes and a degree of integration. The moving average process models the observed time series Zt as a linear combination of a sequence of uncorrelated random variables. The autoregressive process models the present value of the time series as a function of its previous values plus a random effect. Integration is used for time series that do not have a fixed mean to make analyses more reliable. Therefore, if a time series is stationary after d times of first-order differencing and is modeled as ARMA (p, q), then the original time series is an autoregressive integrated moving average time series in the form of ARIMA (p, d, q).


## 3. Results and Discussion

After implementing the Box–Jenkins method in EViews software for the demand time series of one customer, the model ARIMA (0,0,1) was chosen for fitting and forecasting. Therefore, the fit equation for the independent variable demand was obtained as follows:(25)1−B1demt=q0+q1Bat.Using the software, the forecasts for the first difference were obtained as shown in Figure1:Figure 1
Forecasts for the first difference of customer 6.The red, green, and blue lines in this figure indicate the actual, forecasted, and residual values, respectively. The demand forecasts obtained for the first, second, and third periods of the planning horizon are all 37 units.The model presented in Section3 was solved with GAMS software using the complex solver on a computer system with a 2.53-GHz CPU and 8 GB of RAM. The studied supply chain has two levels, with one central warehouse on the first level and 10 customers on the second level. The transport fleet consists of two types of vehicles (8-ton and 10-ton trucks). In the solution obtained after 1000 seconds, the objective function value is 3324.994 (thousand tomans). The optimal route is constructed such that the expected demand is met.Regarding the inventory levels of nodes, only node 6 had kept 1 unit of inventory at the end of period 2. The amount of product received by each node was also determined based on the change in the vehicle load and after visiting the node.In the following, we use the solutions of the model to provide some analyses and interpretations that might benefit managers and researchers interested in the field. If we assume that the travel time does not depend on the time of day in which the trip takes place and also does not define any time constraint for the service, the order and sequence of service will be as follows (Figures2–4):Figure 2
Details of the routes formed in the first period.Figure 3
Details of the routes formed in the second period.Figure 4
Details of the routes formed in the third period.Under this assumption, the total distance traveled will be 375.4 (km), which is less than the 450.75 (km) obtained in the original solution. However, since time-dependent variable costs make up a major portion of the total cost, one can surmise that a reduction in the routing cost can offset the increase in distance traveled. Likewise, considering the social and environmental impacts of decreasing the travel length of distribution vehicles, such as reduced traffic load and thus reduced air pollution and travel time of citizens, managers of distribution companies might be able to make a profit by spending less on the penalties imposed by the government or earning more from the rewards offered by the government for traffic control purposes.As the solutions of the model demonstrate, in certain instances, the model has chosen the longer path between the two nodes. For example, in the first period, the model has chosen the slightly longer second route between nodes 1 and 11 (the 16.5 km long second route rather than the 15.1 km long first route).While this may not conform to the classical view of cost minimization, with a more holistic view of the matter, including the environmental and social benefits that lie within this choice, one can understand the importance of minimizing the total time of tours. This travel time reduction makes it possible to serve more customers within the bounds of the delivery horizon and also reduces time-dependent variable costs such as driver wages, which can make up a large portion of the routing costs.Therefore, managers working in the field of distribution should arrange the order in which customers are served such that those positioned in crowded areas with heavy traffic during peak hours are visited only during off-peak hours (e.g., the middle of the day as opposed to the early hours), when they can be accessed more quickly and via less busy roads. Moreover, this is evident in our solutions in the sense that customers who were located in the city center, where traffic is heavy for a few hours around the start and end of office hours (during peak hours), were served in the middle hours of the day, and customers located on the outskirts of the city were served near the beginning or end of the delivery horizon.As can be seen, costs are almost linearly related to the spoilage rate (i.e., the lower the spoilage rate, the lower the costs). For example, reducing the spoilage rate by half will reduce the inventory routing cost by 1.7%, and reducing it to 0.2 of its current rate will reduce the inventory routing costs by 0.6%.


## 4. Conclusion

The pharmaceutical industry is immense, with a broad range of logistical operations. Considering this industry's ongoing technological and competitive developments, pharmaceutical companies must uphold modern supply chain management principles to maintain their competitive advantages and survive. As in most other supply chains, inventory holding and distribution costs account for a large portion of pharmaceutical supply chain costs. Thus, any cost reduction in these two areas can lead to significant cost savings for these companies. This study modeled an inventory-routing problem for perishable pharmaceutical products with time-dependent travel times and uncertain demand on a multi-graph. The reason for adopting the multigraph approach and defining the travel time as a function of the time of day on which the trip occurs (i.e., traffic volume) was to make the proposed model more aligned with real-world pharmaceutical product distribution networks. The model was formulated to guarantee the FIFO property with a conversion of the travel speed function to travel time. The Box–Jenkins forecasting methodology was employed to control the uncertainty in demand. Given the fact that most drugs have a fixed shelf life, perishability was modeled by determining a variable representing the spoilage rate. To validate the proposed model and demonstrate its applicability, it was applied to the actual data of Asia Pharmaceutical Distribution Company and, mainly, the demand data it has recorded since 2017. The solution results indicated the possibility of achieving cost-saving and reducing product spoilage. The total cost of the supply chain for this particular case was computed to 3324,994 monetary units. Furthermore, adopting this model's solutions brings about some inherent social and environmental benefits, such as relieving traffic load and emission.Regarding future studies in the respective arena, the recommendations below are provided:(1)
Employing environmental and social criteria to analyze the effect of utilizing this type of inventory routing problem on the total traffic load, emission, and product spoilage.(2)
Modeling perishability as the residual life of products based on their age.(3)
Considering other causes of the time dependence of travel time, including weather circumstances and traffic accidents in the decisions and examining their effect on the solutions.



---

*Source: 4526641-2022-09-29.xml*