Supply chain planning of vaccine and pharmaceutical clusters under uncertainty: The case of COVID-19

As an abrupt epidemic occurs, healthcare systems are shocked by the surge in the number of susceptible patients' demands, and decision-makers mostly rely on their frame of reference for urgent decision-making. Many reports have declared the COVID-19 impediments to trading and global economic growth. This study aims to provide a mathematical model to support pharmaceutical supply chain planning during the COVID-19 epidemic. Additionally, it aims to offer new insights into hospital supply chain problems by unifying cold and non-cold chains and considering a wide range of pharmaceuticals and vaccines. This approach is unprecedented and includes an analysis of various pharmaceutical features such as temperature, shelf life, priority, and clustering. To propose a model for planning the pharmaceutical supply chains, a mixed-integer linear programming (MILP) model is used for a four-echelon supply chain design. This model aims to minimize the costs involved in the pharmaceutical supply chain by maintaining an acceptable service level. Also, this paper considers uncertainty as an intrinsic part of the problem and addresses it through the wait-and-see method. Furthermore, an unexplored unsupervised learning method in the realm of supply chain planning has been used to cluster the pharmaceuticals and the vaccines and its merits and drawbacks are proposed. A case of Tehran hospitals with real data has been used to show the model's capabilities, as well. Based on the obtained results, the proposed approach is able to reach the optimum service level in the COVID conditions while maintaining a reduced cost. The experiment illustrates that the hospitals' adjacency and emergency orders alleviated the service level significantly. The proposed MILP model has proven to be efficient in providing a practical intuition for decision-makers. The clustering technique reduced the size of the problem and the time required to solve the model considerably.


Introduction
As epidemic outbreaks are highly life-threatening and require an immediate and proper response, and also considering the growth of disastrous events, it is highly recommended to regard risks as vital elements of supply chains and reserve an emergency supply chain design for this sake [1,2]. After an epidemic occurrence, systems come to a halt, and if emergency operations are not triggered, organizations will face long-lasting suffering. Today, with the spread of the coronavirus (COVID-19), many decision-makers are perplexed about managing their systems. The economy has been highly affected by the consequence of this epidemic, and an economic recession was predicted by the World Economic Forum (WEF), World Bank, and International Monetary Fund [3]. Companies must plan how to recover from the damages that occurred to their systems, especially when knowing that it takes time to recover from this harm.
As reported by the Center for Disease Control and Prevention (CDC) and World Health Organization (WHO), since the start of the pandemic, several variants of the SARS-CoV-2 virus have emerged, some of which are more transmissible and potentially more severe than the original strain. Also, it is known that while COVID-19 vaccines have been shown to be highly effective at preventing severe illness and death, breakthrough infections can still occur. However, vaccinated individuals who get infected with the virus tend to have milder symptoms and a lower risk of hospitalization and death compared to unvaccinated individuals.
The damages striking the healthcare systems due to epidemics are more severe than other organizations because it affects healthcare systems' functionality and changes overall patient demand dramatically [4,5]. To that aim [6], helped to identify health system units and show the capabilities of the resilient unit during the COVID-19 epidemic. One of the most vulnerable sectors of healthcare systems is its distribution network which relies on third parties to perform. Pharmaceutical distribution in healthcare is carried out through the traditional supply chain and cold chain. A cold chain is referred to as a state of preserving perishable products through entire production processes for consumption in a specific range of temperatures [7]. Temperature-sensitive medicines are an integral part of contemporary supply chains; thus, suitable design and control of cold chains are compulsory to reduce waste in monetary and product terms.
With the ever-growing complexity of healthcare supply chains and increasing demand for pharmaceutical products due to the contagious virus spread, more practical approaches are required. Moreover, hospitals need to acquire a supply chain that is responsive during disruptions and can consider pharmaceuticals as a whole, taking into account their different characteristics (e.g., preservation conditions, shelf life, volume, and priority). In this regard, we propose an integrated cold and non-cold multi-commodity supply chain optimization model that involves different pharmaceutical groups, called Pharmaceutical Clusters (PhCs), with unique features representative of constituent pharmaceuticals. In order to analyze the real-world challenges faced by hospitals to procure different medicaments and due to the vitality of this issue, this paper provides the following contributions to the problem.
• A grouping methodology reduces the model size and involves various characteristics of pharmaceuticals. This paper uses the K-means clustering algorithm to consider the high number of various pharmaceuticals and their features; • A multi-commodity mathematical model that integrates features of cold chain and traditional supply chain and deals with the basic nature of epidemics and uncertainty. Integrating the cold chain and traditional supply chain allows for taking into account pharmaceuticals with different temperature preservation conditions. Considering the cold chain and uncertainty both improve the model's applicability; • Addressing the system's cost minimization and service level simultaneously. Assessing these two integral indexes best describes the responsiveness and efficiency of the proposed model; • Considering characteristics mostly neglected in the literature, such as shelf life, and proposing some insights to keep the model responsive during different stages of the COVID epidemic.
The remainder of this study is organized as follows. Section 2 presents a literature review; the problem statement is provided in Section 3; Section 4 proposes the numerical experiment and analyses results induced from the mathematical model, and Section 5 entails conclusions drawn from the study.

Literature review
Analysis of data from Web of Science indicates a notable increase in the number of research studies published in the fields of Healthcare Supply Chain and Cold Chain in recent years, as demonstrated by the upward trend illustrated in Fig. 1. The increasing number of research studies in Healthcare Supply Chain and Cold Chain highlights their growing importance in the healthcare industry. However, more attention and effort are needed to address the challenges and fully realize their potential. Researchers can drive progress by dedicating more resources to these areas, ultimately contributing to better health outcomes and more efficient healthcare systems.
According to the data presented in Fig. 2, there has been a noticeable increase in the size of the cold chain logistics industry, which is projected to continue on an upward trajectory until 2028. This growth pattern emphasizes the cold chain importance and the rising need for more studies in this field.
Dixit & Routroy [8] indicate that supply chain operations, performance measurement, inventory management, lean and agile operation, and application of information technology are thoroughly researched, but employee and customer training, tracking and visibility of medicines, cold chain management, human resource practices, risk management, and waste management demand more extensive investigations. To thoroughly analyze the related former studies, a review has been carried out in the healthcare supply chain, cold chain, and supply chain during epidemic outbreaks.
So far, numerous methods have been used to model different supply chains, and the most popular ones are mathematical modeling and simulation [2,9]. A study carried out by Kelle et al. [10] illustrates a framework to enhance the current pharmacy inventory management policy. This research proposes a mathematical model that uses min and max par levels, resulting in pharmaceutical inventory cost reduction by 70%-80%. In the context of pharmaceutical logistics study, drug variety (formulary) is relatively central.
Sheikh-Zadeh et al. [11] proposed a performance-based inventory classification solution for a multi-item two-echelon (r, q) replenishment model, which proved to be outstanding compared to traditional ways of classification with only a slight error compared to the complete enumeration solution. In addition to inventory management products, categorization can assist supply planning. Kees et al. [12] presented an optimization model that considers three different physical, economic, and information flows among supply chain nodes and provides a multiperiod approach to minimize two interacting objectives of the total cost and total product shortage. They used the traditional ABC classification of inventory items in their MILP model.
On top of the studies mentioned above, there are many pharmaceuticals that are categorized as perishables and need conservation in a specific temperature range. Insulin, some eye drops, some ear drops, Vaccines, Biologicals, and most HIV test kits are examples of temperature-sensitive pharmaceuticals. Unfortunately, these products have been widely neglected in the literature, which makes the existing approaches unpragmatic. In an attempt to analyze shelf life, Priyan & Mala [13] provided a production-inventory game theory model to include a continuous quality degradation rate in a two-echelon hospital and pharmaceutical company supply chain. In a recent study conducted by Shafiee et al. [14] pertaining to the evaluation of risk in supply chain networks involving perishable goods, it was revealed that the perishability of the products in question plays a pivotal role in instigating and exacerbating other risk factors.
Some studies considered environmental issues of the cold chain as it is considered an energy-intensive supply chain, due to refrigeration, with an estimate of using 15% of the energy consumption of the world [15]. Meneghetti et al. [16] modeled a multi-period frozen food transport problem. In their model, they considered traffic congestion, multi-drop deliveries, and outdoor temperature, which contributed to real-world conditions while minimizing fuel consumption for both transportation and cold preservation of products. Mosallanezhad et al. [17] developed a supply chain network for fresh seafood that considers sustainability issues. They proposed a multi-objective mathematical programming model to take care of monetary and waste recycling aspects at the same time. Their results show that the multi-objective golden eagle optimizer algorithm is the best technique to address the computationally challenging model proposed by the authors.
Moreover, Gallo et al. [18] designed a sustainable cold chain for the distribution of food, considering multiple routes and multiple transportation modes. By considering food shelf life and optimum temperature set-point for food preservation, they developed a mixed-integer linear programming (MILP) model to minimize total energy consumption. It is relevant to highlight that Abbasi et al. [19] presented a framework that focuses on achieving sustainable indicators for the supply chain during the COVID-19 pandemic. Tirkolaee et al. [20] developed a bi-objective socio-economic optimization model to tackle the challenges of high costs and job scarcity in the blood supply chain network during the COVID-19 pandemic. The model took into account uncertain factors such as demand, capacity, and blood disposal rates. The researchers demonstrated the effectiveness of their approach by applying it to a real-world case study in Iran.
Another context in which the cold chain gains great value is vaccine procurement, especially during epidemics and outbreaks. The analysis of the supply chain during an epidemic and the damages that epidemic conditions deliver to supply chains are still novel areas of study and one demanding thorough investigation [21]. Now, with the coronavirus (COVID-19) outburst, this demand calls for greater attention from researchers to establish guidelines for managers and decision-makers, especially in the healthcare supply chain field, to minimize the abundance of uncertainties and risks and usher managers to make informed decisions.
Tavana et al. [22] introduced a MILP formulation for ensuring an equitable distribution of COVID-19 vaccines in developing nations. Their approach took into account the temperature requirements of the vaccines, categorizing them into cold, very cold, and ultra-cold types, thereby providing a comprehensive framework for effective vaccine allocation. Kamran et al. [23] conducted a study aimed at analyzing the supply chain network for the COVID-19 vaccine during the ongoing epidemic. The study focused on investigating various aspects of the vaccine's production, distribution, location, allocation, and inventory control within a four-echelon supply chain network. To address this complex problem, the researchers utilized a simulation-optimization approach to develop their model.
To address the various disruptions in the supply chain caused by the rapid growth of the disease, Dasaklis et al. [24] formulated a smallpox outbreak through two modules of a mathematical model. The first module casts the spread of the disease, and the second considers the distribution of multiple medical supplies. It is shown for this particular disease that imminent response (vaccination within four and/or nine days) is required to eradicate the disease because, after 36 days, it is almost impossible to deal with the rapid spread of it. Ebrahimi et al. [25] conducted research on the value of open data systems in coping with uncertainties in healthcare systems during the COVID-19 pandemic. The findings indicated and ranked the barriers to using open data by using interpretive structural modeling and the opinion of eight healthcare experts.
Infrastructures for implementing an efficient vaccine supply are of great value, and different countries face different challenges in designing the vaccine supply network. Designing a cold supply chain network involves considering various factors such as the type of vaccines or medical products being transported, the temperature requirements for each product, the available transportation modes, the infrastructure and facilities available for storage and distribution, and the geographic and logistical challenges of the region. Developing countries often face unique challenges in implementing an effective cold supply chain network, including limited resources and infrastructure, inadequate transportation systems, and difficulty reaching remote or underserved areas. According to Ref. [26]; it is imperative to design a cold supply chain network that takes into account these challenges and incorporates appropriate solutions, such as alternative transportation modes and innovative storage and distribution methods, it is possible to increase access to life-saving vaccines and medical products in developing countries, thereby improving public health outcomes.
Fattahi et al. [27] provided a stochastic programming model to allocate resources to the patients optimally in COVID-19 epidemic conditions. A study was conducted on resource planning strategies that could help health systems effectively manage a surge in demand during epidemics. The strategies that were studied included integrated healthcare resource allocation and sharing, as well as patient transfers. Khodaee et al. [28] addressed the equity problem of COVID-19 vaccine distribution in a cold supply chain network. They used a MILP formulation to model the supply chain planning and applied their work in European countries. In a related context, Mosallanezhad et al. [29] undertook a study on the supply chain of personal protective equipment during the COVID-19 pandemic. To address the issue of cost and unfulfilled demand within the system, they proposed a multi-objective model.
Given that vaccine supply chain risks during COVID-19 can lead to mortalities, Valizadeh et al. [30] proposed a bi-level optimization model that considers both the risks associated with the vaccine supply chain and the monetary aspects of the system. The upper level of the model addresses the vaccine supply chain risks, while the lower level focuses on financial considerations. Moreover, the model accounts for the uncertain demand for vaccines during the pandemic. Sadjadi et al. [31] designed a cold vaccine supply network for the vaccine industry in Iran and considered different uncertainties in vaccine provision and distribution. A MILP model was proposed for the formulation, and a robust counterpart approach was utilized for addressing uncertainties. Yang et al. [32] designed a cold vaccine distribution network for low and middle-income countries, taking into account WHO guidelines for children's universal vaccination. They provided a mixed-integer program (MIP) model to manifest different network attributes and organized a novel algorithm to solve large problems. The model was tested for four African countries.

Research gaps
In order to identify research gaps, Gupta et al. [33] conducted a review of existing research on pandemics and epidemics in the fields of operations management, operations research, and management science. Based on their analysis, challenges were identified, and future research streams were recommended. The study highlights the importance of examining the management of specific diseases within the context of healthcare supply chains. It notes that understanding the unique characteristics of different diseases can help healthcare organizations tailor their supply chain strategies to address the specific needs of patients and healthcare providers.
Furthermore, the study identifies a significant research gap in the area of healthcare supply chains, particularly with regard to hospitals and emergency departments. The authors note that while much research has been conducted on healthcare supply chains in general, there is a lack of attention to the unique challenges and opportunities posed by hospitals. By identifying this gap, the study also delineates that there is still room for contributions in the development of vaccine distribution logistics.
Upon reviewing various research studies, certain key factors have emerged that are crucial in identifying gaps in the existing literature. Some of the most significant of these factors include.
• Considering the cold chain design in developing countries [26]. The study of the cold chain system has the potential to enhance the availability and distribution of vaccines, particularly in developing nations.
• Shelf-life and expiration dates of commodities [13,34]. The significance of ensuring high-quality commodities during distribution has become increasingly paramount. • Integration of planning pharmaceutical items as a whole. The majority of cold chain studies have primarily focused on vaccines as the sole commodity of interest [34]. • Impractical methods for grouping multi-commodity supply chains are used, or supply chains are limited to certain commodities to avoid the complexity of dealing with a large-scale supply chain model. • Cold chain management has received insufficient attention in the existing literature, indicating a significant research gap in this area [35].
Furthermore, Table 1 brings out differences amongst existing studies and this paper and indicates research gaps.

Problem statement
The problem comprises the distribution of medicine groups provided by suppliers within a three-month time horizon in the COVID-19 epidemic condition. The cold chain is also integrated with the proposed supply chain problem to alleviate degradation effects on temperature-sensitive pharmaceuticals. The supply chain network considered in this research effort consists of four echelons, including suppliers, a Demand Hub (DH), hospitals' pharmacies, and hospitals' care units, as shown in Fig. 3. The hospitals are located in proximity to each other in a city, and through a shared real-time hospital information system (HIS), a third-party logistics provider (3 PL) controls the level of inventory in the pharmacies. Almost all of the decisions of inventory control and management are made by the 3 PL. In this way, task division is done properly, and inventory management is assigned to the professionals. This means hospitals focus exclusively on delivering highquality care, resulting in the highest possible service level.
Suppliers are responsible for procuring different medicaments and sending the requested number of drugs to the DH. They have a limited capacity for procurement and can only deliver certain medicine groups. The DH is governed by the 3 PL and is located closer to the hospitals. Due to epidemic conditions, the high demand rate of the whole country for medicines, and the long distance between suppliers and the DH, the 3 PL ordering frequency is limited. Hospital pharmacies are able to store medicines to address the demands of both outpatients and inpatients. Additional drugs can be delivered within a lead time by the DH. Moreover, unfulfilled demand for care units can be settled through emergency orders-supplied from other pharmacies in the hospital's vicinity. Ultimately, the main goal of the proposed study is the minimization of the underlying system's cost.
The basis of the network design is inspired by the Supply-Demand Hub in Industrial Clusters (SDHIC) model proposed by Kayvanfar et al. [47] and Kayvanfar et al. [58]. The supply chain network considered in the aforementioned study is proved to be the most efficient compared to other proposed models; However, this network is best suited for industrial supply chain networks, including manufacturing. The network proposed in our model encompasses the characteristics modified explicitly for pharmaceutical supply chains inspired by Iannone et al. [48]. Iannone et al. [48] proposed an integrated supply chain for drug distribution consisting of a central pharmacy, hospital pharmacy, and medical unit. The information about patients' prescriptions was deemed to be shared through electronic medical records (EMR) with other network participants. The advantages of the network designed in the present study could be presented as follows.
• Inventory cost mitigation as 3 PL stores medicines in a centralized warehouse and manages the orders; • Reduction of inventory specialists and practitioners at the hospitals; • Increasing hospitals' efficiency as they only focus on delivering medical services, which during epidemic conditions can reduce work pressure; • Capital investment reduction as pharmacies store a fewer number of medicines; • Cost reduction of pharmacies pertaining to lower holding and labor costs; • Risk reduction, owing to emergency deliveries from other pharmacies.
It is worth mentioning that the name of the "Demand Hub" is inspired by the network, as it is closer to the hospitals and helps manage and control the demands of the patients effectively.

Assumptions
The following assumptions are made in the present study.
• Each time frame in set T equals one day.
• Set T has 90 elements to demonstrate three months' horizon planning of winter as demand changes due to epidemic volatility in each season.
Considered supply chain network, including suppliers, DH, pharmacies, and care units.
• The lead time is deterministic and equals one day for every order except emergency orders, which are negligible, and their lead time is deemed zero due to the proximity of hospitals. • The expiration rate for each pharmaceutical group is known, and excessive storage causes degradation over time. • The demand for each cluster is deterministic throughout the time horizon and is achieved using historical data. This assumption is novel as it is inspired by hospital supply chain practitioners that use seasonal historical data to estimate the demand parameter. • Some pharmaceuticals, like the COVID-19 vaccine, are deemed to be apportioned, and their procurement is not based on demand. This assumption is novel and congruent with the early stages of vaccine distribution when vaccines are not largely accessible. • Inventory management of pharmacies and DH is dedicated to the 3 PL. As mentioned by Kayvanfar et al. [47], hospitals, such as SMEs, can experience hardships relating to profitability, business planning, low efficiency, decision-making information, etc. • The supply capacity is considered to be limited due to epidemic conditions. • Priority is included in the proposed model as PhCs with higher priorities have higher shortage costs to prevent vital pharmaceutical deficiency. This assumption is almost unprecedented in supply chain studies. • Each supplier can only provide certain groups of pharmaceuticals, and orders are settled once a month. • The shortage is considered in the model and is lost.
• A higher transportation and expiration cost is considered for cold pharmaceuticals. This novel assumption incorporates the high cost of preservation and distribution of commodities in the cold chain. • Trucks with different capacities are considered in the model as means of transportation. The refrigerator trucks are included to deliver temperature-sensitive pharmaceuticals. • Demands of every pharmaceutical in a cluster are averaged to reduce the model complexity and provide decision-makers with a flexible and nimble model to test different assumptions during the epidemic when time-consuming models are impractical, and swift decisions are required. This assumption is unexplored in supply chain studies.

Model formulation
The definition of notations included in the following mathematical model is presented below. It is obvious that the number of sets is reduced due to the use of the clustering approach before modeling, consequently resulting in the model's size reduction. The impact of implementing this approach will be further discussed in the next section. The approach used for clustering the pharmaceuticals is an unsupervised learning method called K-means clustering. The algorithm is completely explained in the consecutive sections, but it is worth noting that this approach assigns pharmaceuticals into clusters based on their proximity of pharmaceuticals' features, including temperature, demand, priority, and volume. The outcome of the clustering method shapes the PhCs defined in the mathematical formulation. The proposed MILP model is given as follows. The objective function of this model is to minimize the overall supply and distribution costs pertaining to the hospital supply chain.
The contributions from the first term and second term of the objective function (1) account for total storage cost at the DH and the pharmacies, respectively, through the time horizon. The third term attributes to the shortage cost. The successive two terms pertain to the expiration costs of pharmaceuticals at the DH and the pharmacies. The pharmaceuticals' purchasing costs to the DH and the pharmacies are associated with the sixth and the seventh term. The eighth term illustrates the emergency purchasing cost of care units from other pharmacies. The last three terms manifest the transportation costs from the suppliers to the DH, the DH to the pharmacies, and the emergency cost of transportation from the pharmacies to the other care units.
The constraints of the model are presented as follows: Constraint (2) balances the storage level at the DH for each PhC at period t. It ensures that the number of pharmaceuticals stored is not expired and indicates that the inventory level equals the inventory level from the last period plus orders received from the suppliers minus orders sent to the pharmacies.
Constraint (3) is also an inventory balance constraint, except that this one is for the pharmacies. It makes sure that the stock amount at each period equates to that of in the last period, adding the incoming flows from the DH and subtracting the outgoing flows to the care units. ∑ Constraint (4) indicates that pharmacies deliver PhCs in each period based on the demand of patients. Provided that the demand of patients at each period may not be satisfied, the unsatisfied portion is accounted as a shortage. By assigning a high cost for shortage, the model prevents shortage accretion. ∑ m X k,m,i,t ≤ Ca k,i ∀k, i, t Constraint (5) assures that at each period, the quantity of PhCs delivered from each supplier is less than the capacity of that supplier. ∑ m X k,m,i,t * Vol k ≤ SCb k * r k,i,t ∀k, i, t Constraint (6) regulates the deliveries to the DH as supplies from suppliers are allowed as long as the suppliers hold the ability to procure a specific PhC. Furthermore, this constraint indicates that the maximum procurement volume from each supplier at period t equals the volume capacity of the DH for that particular PhC.
Constraint (7) determines the ability of each supplier to procure each type of PhC at each period. In other words, the supplier i can procure the PhC-k at each period if the parameter ps k,i equals 1.
Constraint (8) specifies the frequency of placing orders by the DH from each supplier and for each PhC. The DH is allowed to place orders once a month.
Constraint (9) enforces that the volume of each PhC received at each pharmacy should be at most equal to the storage capacity of that pharmacy at each period.
Constraint (10), just as constraint (8), delineates the order frequency of each pharmacy to be equal to ten days for each PhC.
Constraint (11) makes sure that the outgoing flows to the pharmacies from the central hub are at most equal to the stored PhCs inside the DH plus those received from the suppliers. In better words, this constraint balances the available pharmaceuticals in the DH with the outgoing pharmaceuticals to the pharmacies. ∑ Constraint (12) indicates that the maximum quantity of pharmaceuticals delivered to the care units equals the number of pharmaceuticals available on hand. Considering this constraint with the last one makes it clear that the balance throughout the whole supply chain is met.
Constraint (13) expresses that the volume of PhCs at each period in the DH cannot exceed the size of the DH storage.
Like constraint (13), constraint (14) assigns a limit for the volume of the stored PhCs at the pharmacies. It indicates that, on aggregate, the volume of pharmaceuticals stored at each pharmacy should be less than the total capacity of that pharmacy for that particular PhC.
Constraint (15) determines the number of vehicles required for transporting PhCs from suppliers to the DH at each period considering the capacity of each vehicle. This constraint specifies that the number of vehicles required is at least equal to the accumulative volume of the pharmaceuticals supplied over the capacity of that particular vehicle type.
Constraint (16), just as constraint (15), makes known the number of vehicles needed to transport deliveries from the DH to the pharmacies.

Numerical experiment
In this section, a case study of three hospitals in Tehran is considered to validate the proposed MILP model and the proposed approach's capabilities. For this sake, section 4.1 presents data collection, K-means Clustering, and model characteristics. Section 4.2 analyses the flow of the supply chain by giving an illustrative example. Section 4.3 reports the results drawn from the case study. Section 4.4 presents the stochastic solution of the model. Section 4.5 depicts the sensitivity of particular variables to changes in some parameters, and section 4.6 offers some managerial insight based on results drawn from the case study.

Data collection, K-means clustering, and model characteristics
The numerical experiment carried out is based on clusters drawn from the clustering method. The case study includes three hospitals in Iran's capital city, Tehran, and the assumptions are fit for this case. The data is extracted from the Iranian Food and Drug Administration. The dataset used in this study is drawn from the Annual Pharmaceutical Statistics of March 2019 to March 2020. Additionally, the statistics were filtered to the winter of 2020 to account for this season's data concurring with the first surge in the number of confirmed cases of patients with the COVID-19 virus.
In order to control a large number of commodities, researchers have proposed many techniques that help group commodities. One of the most popular methods used for this purpose is the ABC analysis [12,49]. This classic technique follows the Pareto 80/20 rule and groups commodities so that higher-value commodities with lower volume are higher in rank. The developed techniques have been initially used for inventory control and management studies, but due to their practical application, they can be used in mathematical modeling as well. Additionally, each medicine can have various features like size, volume, priority, storage condition, etc., that may perplex decision-makers to analyze them, and classic techniques fail to consider them.
Methods using K-means clustering algorithms are more competitive than the traditional ones [50], yet they are mostly ignored in practice; as to the extent of our research, no study utilized this technique to group commodities for the mathematical model. K-means clustering is a technique developed by MacQueen [51] that clusters data points based on their affinity. It has been used by numerous recent studies for different applications [52,53]. The advantage of this technique over traditional techniques like the ABC analysis can be.
• All features in a multi-feature data can be included in K-means Clustering; • The K-means clustering technique is less subjective and less dependent on expert opinion, and clusters are shaped only based on the affinity of items to each other; and • K-means clustering may need more computation but is implemented significantly faster for big data.
The K-means clustering algorithm is a machine-learning technique that identifies groups of data points in a multidimensional space [54]. Through the implementation of this algorithm, the number of groups, as well as the groups themselves, are determined. What K-means does is to take out the complexity from the grouping process and help identify clusters of data points in our dataset through an iterative algorithm. The intuition of this algorithm is given in Fig. 4.
The algorithm starts with selecting the number of clusters. For this purpose, a metric is used, which is called Within Cluster Some of Squares (WCSS). This metric calculates the optimal number of clusters for a dataset by solving different K-means clustering problems with different values for the number of clusters (K). For each clustering problem sum of squared distances of data points in a cluster from their central point, called the centroid, is calculated. By increasing the value of K and converging to the number of data points, the algorithm assigns each data point a centroid on itself that equates WCSS to zero. However, after a certain number of iterations, WCSS does not change significantly, which means that from that point on, increasing the number of clusters is computationally redundant.
Therefore, different WCSSs are plotted on a chart against the different number of clusters, and through using the Elbow Method, the optimal number of clusters is identified. After finding an optimal K, a random point for each cluster is selected as a centroid. The algorithm then continues until centroids are assigned to a point in space where their distance to their cluster members is minimized, and no readjustment for centroids can be done. Numerous methods can be utilized to calculate distances, but the most famous method is Euclidean distance, which is also used in this study.

K-means algorithm implementation
The data set used in this research is gathered from the pharmacies of three hospitals' daily sales reports. The data set comprises 562 unique medicines and integrates the daily demand of these pharmacies over a three-month period in the winter of 2020. The data has four features quantity, priority, volume, and storage temperature. Storage temperature required for maintaining the medicines' potency falls into two groups Controlled Ambient Storage (CAS)-from +15 to +25 • C-and Chilled Storage (CS)-from +2 to +8 • C.
The data is first averaged over daily demand and categorized based on the aforementioned groups of storage temperature (as medicines in both groups cannot be clustered together due to degradation). Then each category is clustered using the K-means algorithm. From over 562 different medicaments requested from different pharmacies, only 45 of them demand CS. This category does not get clustered because the pharmaceuticals requiring CS are low in number and need to be delivered together at the same temperature. On the other hand, medicines stored in CAS are clustered in three different groups according to Fig. 5the Elbow method analysis. In order to avoid the random initialization trap and make faster convergence, the K-means++ extension developed by Arthur & Vassilvitskii [55] was executed. The result is reported in Table-2. Overall, the purpose of this study to cluster medicines is threefold.
1. To have better control on medicines, as 562 individual types of medicine are difficult to track and control.  It is worth mentioning that all pharmaceuticals grouped in a cluster are treated fairly, and all the processes, including supply, storage, control, and transportation, happen the same way. Also, the model parameters are chosen according to the features of these clusters. For instance, the PhC with higher priority has a higher shortage cost. This higher cost acts as a penalty for the objective function and prevents the system from not satisfying high-priority demands (see Table 3).
As mentioned previously, the implementation of the clustering algorithm substantially reduces the number of constraints that affects the speed of the solver significantly. To drive the point home and illustrate the functionality of this approach, it is required to consider the number of constraints without clustering. The comparison of the differences in the sets is drawn in Table 4. As it is evident from the table, if pharmaceutical clustering is eliminated, some sets need to be appended to the model, and an increase in the members of the pharmaceutical set is inevitable. By taking the changes of sets into account and using a rule of thumb to calculate the number of constraints, it is concluded that the number of constraints at least rises to 44,904,878 individual constraints.
Knowing that constraints are a function of O(h 2 ) and have a proportionally high effect on the model-solving time, the clustering approach has proven extremely useful in multi-commodity multi-characteristic problems. Although it may seem that the original problem may distinguish from the problem with clustered pharmaceuticals, it is considered a cogent trade-off, knowing that the original problem is practically impossible to solve. Given the epidemic conditions, the clustered model will satisfy the policy-makers' urges to reach an overall understanding of the system.

Model parameters
The studied supply chain involves three hospitals with three local pharmacies inside the hospitals (c = 3, d = 3). The pharmacies are procured by a 3 PL who is in charge of managing the DH as well. Overall, 562 pharmaceuticals are delivered to the DH in four PhCs (k = 4) through four suppliers (i = 4). PhC-A is available through Supplier-1 and Supplier-2, PhC-C is delivered by Supplier-1 and Supplier-3, and PhC-B and PhC-D are only procured by Supplier-2 and Supplier-4, respectively. A planning horizon of three months (T = 90) is considered, which corresponds with three months of winter. Moreover, it is assumed that the pharmaceuticals requiring CS are expired on average after 60 days with a specific rate derived from historical data (θ = 0.015). The demand parameter of both outpatients and inpatients is drawn from Table 2. Some of the remaining parameters are drawn from Kayvanfar et al. [57], and the others correspond to the real case. The following table indicates the values of these parameters.
The MILP model is coded in GAMS 25.1.3 and solved via Gurobi on a PC Intel Core i3, 2.53 GHz with 3 GB of RAM. The instance of the model consists of 29,509 single equations and 46,657 variables that encompass 44,137 continuous and 2520 binary variables. The obtained solutions give us an insight into the capabilities of the supply chain to provide necessary medicaments during the epidemic. Furthermore, it provides valuable information on supply chain costs and, more importantly, cold chain costs, enabling decision-makers to make more informed decisions. The solver reached the optimal solution after 2 h and 47 min, which is deemed decent as the model is capable of generating the best solution in a relatively short amount of time for a long planning horizon. This capability enables decision-makers to react to epidemic conditions as soon as possible.

Schematic example
In order to visually demonstrate the flow of pharmaceuticals throughout the supply chain, Fig. 5 presents an illustrative sample of the supply chain. As shown in Fig. 5, the flow of pharmaceuticals requiring   Table 3 The difference in the number of set members with considering clustering or not.   (5,8,10,12) CS and CAS is separated to showcase the cold chain flow. This figure depicts the supply chain with four suppliers, a DH, and three hospitals with local pharmacies, as mentioned before, and the flow of pharmaceuticals from day 60-70. The numbers shown with the black color in this figure are representative of pharmaceuticals with CAS condition, and the ones with the blue color indicate that of CS condition. According to the optimal solution presented in Fig. 6, The model uses both regular and emergency orders to satisfy the demand of the patients. As an instance, the demand of the second hospital from day 60 to day 70 is satisfied by ordering 3630 ambient PhCs and 93 cold PhCs from the local pharmacy and 1898 ambient PhCs and 112 cold PhCs from other pharmacies. Furthermore, the procurement of the DH was carried out on days 61, 63, and 69, which subsequently were added to the stored pharmaceuticals and delivered to pharmacies. All in all, This Figure helps to comprehend better the interactions among the levels of the supply chain.

Numerical results
In this section, the optimal cost of the supply chain is presented, and some further analyses of the optimal variable values are given. The obtained optimal cost by solving the proposed MILP model is calculated to be $2,433,734.951. An upper bound for the model can be estimated to be 4,941,140.92 by leaving all the demands unfulfilled. The cold chain accounts for approximately five percent of it, with a total cost of 121,521.83. Inventory, shortage, purchasing, and transportation costs are included in the objective function, and their optimal values are presented in Table 5. Interestingly, inventory costs account for the majority of the Total Cost.
Just over 32% of the inventory costs pertain to the DH due to the Just-in-Time strategy carried out by the 3 PL, and the rest relates to the pharmacies' inventory levels. Over 13% of the costs are related to the shortage cost, considering the epidemic condition and the abundance of demand. The purchasing cost is the third highest cost pertaining to the supply chain costs. Expiration cost and Transportation cost comprise only a minority of the total cost.
The optimal results for the inventory level of the DH are presented in Fig. 7. It is evident, at first glance, that number of orders did not exceed the limitation posed by equation (8) that there should be only one order every 30 days. The first orders were placed as soon as initial inventory levels were sent to the pharmacies on day 12. Furthermore, PhC-B was the most dominant cluster held in the DH to satisfy the exceeding demand of the patients, yet it did not violate the capacity constraint in equation (13). Fig. 8 illustrates the inventory level of PhCs in the three pharmacies, while Fig. 9 indicates the level of pharmaceuticals requiring CS and are transferred through the cold chain. The limitation of one order per ten days set by equation (10) is not violated, and the level of inventory in Pharmacy C is greater than that of Pharmacy A and B. In addition, it is obvious from Fig. 8 that the level of PhC-D that requires CS is not always the same for the three pharmacies. It can be inferred from this chart that in periods when at least one pharmacy is out of PhC-D, the demand of its care unit is not fulfilled or is procured through emergency purchases from adjacent hospitals.
The number of unfulfilled demands is given in Fig. 10. It is apparent that PhC-B is the primary reason for unsatisfied demands, and tackling its issues would significantly reduce the lost sales. The number of unfulfilled demands has risen over time, starting at around 500 lost sales on day 12 and peaking at just below 1400 pharmaceuticals by the end of the period. PhC-D had a minor effect on total unsatisfied demands, accounting for 64 unsatisfied demands.
In relation to emergency orders placed by the hospital care units, Table 6 presents the number of pharmaceuticals transferred to the care  units from the other hospitals. This alternative way of procurement has enabled hospitals to reduce lost sales and subsequently promote the service level of hospitals. Care Unit B has been benefiting from emergency purchases the most. Overall, around 10,139 PhCs were purchased using emergency orders. Care Unit A and Care Unit C had almost the same number of emergency orders received, with Care Unit A ordering about 600 more PhCs than Care Unit C. It is worth noting that Care Unit A has received the most PhCs with CS conditions throughout the 90 days   period. Table 7, on the other hand, represents the pharmacies capable of procuring emergency orders. In contrast to Care Unit B, which is responsible for putting the highest emergency orders, Pharmacy C is able to provide the most PhCs in case of emergency. Furthermore, Pharmacy B is shown to be the most reliable for transferring pharmaceuticals with CS conditions through the cold chain when other hospitals suffer from a lack of PhC-D.
The number of transportation means used can be observed in Table 8. It is evident from the table that Transportation II had a higher rate of transferring pharmaceuticals, and on aggregate, 355 of them were used to facilitate transportation. In sharp contrast to this, Transportation Mean III has been deemed an expensive mean as it is utilized least during the planning horizon. The rationale for employing diverse transportation modes for pharmaceuticals is evident. There are financial incentives to utilize low-capacity transportation modes for smaller orders, such as emergency requests while employing high-capacity means for larger orders. In addition, as the transportation means are equipped with refrigeration due to the CS and CAS requirements of the pharmaceuticals, selecting transportation means with proper capacity reduces the environmental impact.
In order to study the viability of the freight consolidation strategy, the integrality constraint on transportation means variables was relaxed, and the model was solved again. The number of transportations means used is rounded up, and it is shown in Table 9. This is due to checking the importance of considering freight consolidation in the supply chain. The results after relaxation indicate that a capacity equivalent of over 56 vehicles of each transportation means used was vacant throughout the planning horizon. Although this number seems proportionally considerable, it does not significantly affect the total cost, and it makes the transfer faster and causes diminished lead time. As the cost of unsatisfied demand is higher than using more transportation means, it is reasonable not considering freight consolidation (see Table 10).

Stochastic solution
Uncertainty is an integral part of epidemics. The system's parameters are highly volatile during epidemics, and a minutiae change in the system's parameters is far enough to alter the system's assumptions and invalidate the developed models. In supply chain studies, the most sensitive parameter that reflects the changes in the conditions is the demand parameter. As the epidemic conditions get more severe and the rate of diagnosed or hospitalized patients increases, more and more pharmaceuticals are required to maintain a high service level, and emotional and nominal demands rise significantly.
The uncertainty in parameters arises from experiencing a unique epidemic condition. Hence, historical data used by other researchers are not a valid source to account for the newly arisen outbreak. In order to address uncertainties in our model, before vaccination is widely available for the population, a Wait-and-See (WS) approach has been used. Uncertainty is prevalent before vaccinations are widely available because numerous unknown variables are existent. Therefore, as no historical data can be found, models can consider uncertainty based on worst-case scenarios or wait until the distribution of the uncertain parameters is determined.
The WS approach is viable in our case as historical data on two consecutive winters before the vaccination became widely available in Tehran is accessible. The first winter (the winter of 2020) coincided with the beginning of the first surge in the number of patients diagnosed with COVID-19 in Tehran and attributed to the first scenario of the model. The second winter (i.e., winter of 2021) coincided with the end fourth surge in the number of patients in Tehran and contributed to the second scenario of the model. Based on historical data, the demand parameter in the second scenario is greater than the first by 20%, on average.
The theory pertaining to the WS approach is given as follows [56]: Equation (19) shows the stochastic model with variable x and uncertain variable ξ that is based on scenarios. According to equation (20), in order to achieve WS solution, each scenario's solution must be obtained and calculate the arithmetic weighted mean based on the scenario's probability distribution. The mentioned winters happened before the wide availability of vaccines, and both had the same duration. Hence, their occurrence had the same probability, and the weight of each one of them is calculated to be 1 2 . The optimal solutions of the two scenarios are listed in Table 9.
According to the optimal values in Table 9, the WS solution equals $2,707,604.748. This means that in the long term, the system cost is equal to this value, while vaccination is not widely available.

Sensitivity analysis
An extensive sensitivity analysis is carried out to demonstrate the effects of essential parameters on the pharmacies' total cost and service level. These parameters considered for sensitivity analysis consist of the expiration rate, shortage cost, and storage capacity of the DH. The following subsections pertain to the sensitivity analysis of these parameters individually. Fig. 11a and b represent the effect of changing the expiration rate parameter value on the total cost and shortage value, respectively. It is evident from Fig. 11a that an increase of 0.09 in the expiration rate can impose a $21,690 additional cost on the supply chain. In addition, the rise in the total cost is diminishing as the expiration rate increases. Fig. 11b indicates that the least shortage-or the highest service level-is achieved when the rate equals 1.5%. Any deviation from this rate would cause more shortages. It is worth mentioning that when the rate increases, higher shortages will occur, especially in PhC-D, compared to its decrease.

Analyzing the impact of shortage cost on total cost and shortage value
The impact of changing the shortage cost parameter mirrored in the total cost and shortage value can be seen in Fig. 12a and b. As solving time takes about 3 h, three sets of values for the parameter are considered to draw sensitivity analysis results. In relation to the impact of shortage cost on the total cost, it is obvious that from a certain point, the total cost is indifferent to the changes in the shortage cost. However, when the set of values of shortage cost is decreased by 50%, the total cost experiences an above 32% decrease. Like other trade-offs, this decrease, however, comes at a price. By taking a glance at Fig. 12b, it becomes clear that reducing the shortage cost would dramatically increase the shortage values and consequently causes a considerable decrease in the service level of the pharmacies. In direct contrast, the rise in the shortage cost would maintain almost the same total cost while increasing the Fig. 11a. The change in total cost due to expiration rate variation.   service level to an almost perfect level.

Analyzing the impact of the storage capacity of DH on total cost and shortage value
According to the results drawn from the changes in the value of DH capacity, the status quo of the parameter value in the model is considered the best value, considering the total cost and the shortage numbers. More specifically, according to Fig. 13a, adding or reducing the DH capacity would cause a rise, however negligible, to the total cost. Furthermore, it is apparent in Fig. 13b that changes in the DH capacity would exacerbate the service level.

Managerial insights
Based on the previous sections' model solution results, the following managerial insights can be provided. The aim of providing these insights is that decision-makers and practitioners seeking to implement the proposed model in a similar situation could better understand the context and modify the parameters to fit their conditions the best. Some significant insights are as follows.
1) In order to deal with a proportionally big size of the pharmaceutical supply chain that has intrinsically multi-commodity features, another possible answer to deal with long solving time besides devising competent algorithms is to use grouping techniques that are available to reduce the complexity of the mathematical models. The K-means clustering technique proved beneficial in grouping various pharmaceuticals, and it accounts for the different features of each one of them. Although, it is worth mentioning that these grouping techniques provide solutions that approximate the original problem.
In situations that are as drastic as epidemic conditions, these techniques are worthy of giving practitioners some perspective while making vital decisions in a short time. In situations where the emergency is a non-existent factor, or more precise solutions are required, a higher number of groups can be used on discretion or not used at all. 2) Unlike other conditions where service level does not affect individuals' lives, pharmaceutical delivery during epidemic conditions is of extreme importance. Hence, practitioners should always choose the values that increase the service level for selecting the appropriate parameters' values. For instance, dedicating a high cost for shortages proved to be an efficient way of doing so in the sensitivity analysis section. Moreover, if the condition is not severe and the disease resurgence is not a concern, decision-makers may also consider monetary aspects. 3) In pharmaceutical supply chains where delaying the procurement of medicines could cause irreversible damage to human lives, lateral transshipment that postpones the delivery is deemed disadvantageous. The costs that it would impose on the system are negligible compared to the shortage costs.

Conclusion
Epidemic conditions are among the risks mostly neglected in supply chain policies as they are considered scarce. However, its occurrence widely affects supply chain systems, especially pharmaceutical supply chains. Pharmaceutical supply chains play a critical role during the epidemic as they account for the costs involved and the quality of care. The recent COVID-19 epidemic highly affected these supply chains by skyrocketing the demand for patients and causing an alarming condition for hospitals. In this paper, a decision-making MILP model is provided to aid the procurement planning of hospitals in the COVID-19 epidemic and deal with the characteristics of the various pharmaceuticals, including vaccines. It minimizes the cost involved in a four-echelon pharmaceutical supply chain, namely, ordering cost, holding cost, shortage cost, expiration cost, and transportation costs. The fourechelon supply chain consists of suppliers, a DH, pharmacies, and care units. The DH is managed by a 3 PL responsible for managing the order to and from the DH.
A K-means clustering algorithm was used in this model to group the various pharmaceuticals and vaccines associated with this supply chain. This method proved worthwhile as it enabled the model to consider various features like temperature, priority, and volume and make the model faster to solve, which is essential in an epidemic. It helps reduce the size of the problem by grouping the pharmaceuticals and providing practical implications for practitioners and decision-makers to deal with this highly volatile situation in a proportionally short period. Additionally, based on the current requirements of some pharmaceuticals and vaccines that are temperature sensitive, cold chain characteristics and non-cold ones are considered simultaneously to provide an integrated approach.
Real data from three hospitals in Tehran was used to test the model. The results paved the way for planning an integrated cold and non-cold chain considering numerous pharmaceuticals, including vaccines. The design of the supply chain proved to be responsive in reducing cost and simultaneously maintaining high service quality. Allowance of Fig. 13a. Total cost changes due to DH capacity variation. emergency orders from adjacent hospitals enabled care units to procure the patient's demands optimally. Ignoring lateral transshipment in this research enabled the system to be more responsive during emergency conditions of an epidemic and did not prolong the lead time. The stochastic solution gave an insight into the long-term cost of the supply chain, and it indicated the importance of considering uncertainty. The sensitivity analysis shed light on selecting the correct values for parameters based on hospitals' policies and conditions. If there is a peak in COVID cases, it is recommended to select the values to boost the service level, and if the condition is stable, inserting the values that also consider financial aspects is not prohibited.
The contributions put forward by this study, including integrating pharmaceuticals considering temperature-sensitive pharmaceuticals, integrating cold and non-cold chains, using a grouping technique to tackle the immense size of the pharmaceutical supply chain, and considering the feature of more than 560 pharmaceuticals, are mostly unheard-of in the context of pharmaceutical supply chain and few studies have considered them. Limitations of this study and future research streams are presented below.
• Due to limitations in the available data sources, the present study is limited to hospitals in the city of Tehran. Further studies are required to implement the proposed research in different cities and compare the results. Also, future work can promote this study to a higher level of the outbreak and consider the global pharmaceutical supply chain in a pandemic. • For achieving a maximum level of integration in the hospital supply chain system, the supply chain of other healthcare items, especially fusion pumps, can be considered. This inclusion would further increase the model's applicability. • Euclidean distance was used to measure the distances in the K-means algorithm. Other distances or even other unsupervised learning approaches can be utilized and compared with the current results to evaluate the accuracy. Furthermore, the model can be solved using different unsupervised learning methods, and the solutions can be compared with the actual problem in future studies. • In the uncertain model, it was assumed that scenarios are known, and consequently, the wait-and-see approach was used to achieve the optimal solution. Future studies could consider scenarios unknown and develop methods like Robust Optimization that are not dependent on scenarios and report the results. Also, uncertainty in some other parameters could be considered to capture more indeterministic features of the epidemic.
The configurations employed in this paper and the mathematical model proposed are intended for hospital or pharmacy managers and supply chain practitioners to support their decisions while there is little time to go through the extensive feature of the outbreak. Manipulation of the parameters enables the use of the model during an emergency epidemic or a stable situation.

Data availability
Data will be made available on request. awards, such as ARAP scholarship from A*STAR in Singapore in 2016, and Iran National Elite Foundation during 2012-2016.
Prof. Alireza Haji is a professor of Industrial Engineering at Sharif University of Technology, Tehran, Iran. His main research interests mainly fall within inventory planning and control, project control, stochastic processes, designing the operational systems, and queueing theory. He also has published several papers in peer-reviewed journals and conference proceedings.