A data-driven robust optimization model by cutting hyperplanes on vaccine access uncertainty in COVID-19 vaccine supply chain.

The worldwide COVID-19 pandemic sparked such a wave of concern that made access to vaccines more necessary than before. As the vaccine inaccessibility in developing countries has made pandemic eradication more difficult, this study has presented a mathematical model of a sustainable SC for the COVID-19 vaccine that covers the economic, environmental and social aspects and provides vaccine both domestically and internationally. It has also proposed a robust data-driven model based on a polyhedral uncertainty set to address the unjust worldwide vaccine distribution as an uncertain parameter. It is acceptably robust and is also less conservative than its classical counterparts. For validation, the model has been implemented in a real case in Iran, and the results have shown that it is 21% less conservative than its classical rivals (Box and Polyhedral convex uncertainty sets) in facing the uncertain parameter. As a result, the model proposes the construction of two domestic vaccine production centers, including Pasteur Institute and Razi Institute, and five foreign distributors in Tehran, Isfahan, Ahvaz, Kermanshah, and Bandar Abbas strategically.

Introduction

First reported in Dec. 2019 in Wuhan, Hubei Province, China, and spread rapidly worldwide, COVID-19 is an infectious disease [1] which was declared as a global pandemic by the WHO on March 11, 2020 after affecting 213 countries [2]. The second virus type, much more contagious than the earlier ones, has significantly affected the world industry and economic [3,4] reducing the GDP of industrialized countries and causing the loss of many jobs in in the world [5].

While the disease is spreading pandemically around the world, finding its vaccine is a hard but valuable task that gets harder when millions of people are eager to receive it. The fair delivery and mass vaccination under different constraints requires a resilient robust supply chain network design (SCND). Although the corona pandemic and its double pressure on countries have urged governments to respond better to their citizens by equipping the infrastructures of their health systems, an optimal vaccine SCND that can meet the global needs and provide equitable access for all countries is not an easy task [6].

High inter-company competitions and fear of failure and elimination from the market caused the formation of a new concept - SC management - in the 1980s [7] which emerged the integration of activities because of the boosting number of companies; this inter-company activity coordination gradually led to a revolution in this field [8]. The SC management simply means the integration and coordination of materials, information and financial flow in the chain [9]. Another definition, based on what the Professional Council of Supply Chain Management believes, is that it is the process of planning, implementing and controlling the SC activities in an efficient way [10]. However, the globalization, technology development, companies' fear of failure in the current competitive environment, short life of products, customers’ varied needs/behavior, and most importantly, emergence of diseases with global pandemics and vaccine supply have all forced companies to put the SC management on their agenda and start large investments on it [11], [12], [13].

The SCND, which generally involves making strategic decisions, contains that part of the SC management process that underpins the SC's physical structures and infrastructures. Facility locating has been a widely used SCND aspect in the last two decades; in a review article, authors of [14] divided its related models into continuous, netted, analytical and discrete most of which had commonality in meeting demands and locating facility centers [15].

Optimization of practical problems that requires using real data, always faces parameter uncertainty that can either be stochastic or erroneous due, generally, to the parameter-related lack of knowledge; for instance, demand may be erroneous when determined depending on how precise the problem is physically set in the real world. In such a case, optimization is either stochastic, where the distribution function of uncertain parameters is assumed to be known and all the mathematical modeling relations are based on it [16], [17], [18], or it is robust, where the assumption is not true because the information is so scarce that the decision makers have no idea of the distribution function; here, they generally use a set of realistic values ​​called the “uncertainty set”. Regarding uncertain parameter, basic robust optimization approaches place strict-constraints and do not allow the mathematical model to violate them anywhere in the intended uncertain set. Robust optimization is popular because it performs well in different problems [19], [20], [21].

In the stochastic approach, the key role is played by the distribution function, but in the robust method [22,23], the “uncertainty set” plays the main role in finding more appropriate solutions; the better is the set selected, the better will be the solutions; hence, set selection is vital to the quality of the solutions of the mathematical model. An appropriate “uncertainty set” is that obtained based on various hypotheses and empirical validations.

Recent years’ information explosion and expansion of new data-collection technologies have revolutionized dealing with uncertain parameters. The great access to the information collected from different industries is quite helpful in health systems which are usually affected more due to the uncertain nature of such parameters as the accessibility, demand and system efficiency. Data can make a significant change in robust optimization, which deals with a set of uncertain parameters, and enable, more than before, the implementation of realistic and sustainable cases against the fluctuations of these parameters.

This study that presents a multi-level, multi-period COVID-19 vaccine distribution supply chain MILP model, involves 2 innovations: 1) operational, where the economic, environmental and social dimensions are considered in objective functions 1, 2 and 3, respectively, to create a sustainable model, and the chain addresses both the internal and external supplies, and 2) theoretical that deals with the uncertain accessibility of foreign vaccines. It presents a data-driven robust optimization (DDRO)-Cutting planes approach in agreement with uncertain parameters where the “uncertainty set” is polyhedral to provide enough flexibility for the readers and stakeholders in the intended space.

The rest of the paper has been so organized as to review the literature in Section 2 (focusing on the robust optimization of the vaccine distribution SCND model under pandemic conditions), implement the sustainable SC in Section 3, describe various forms of the proposed approach, using cutting planes, in Section 4 (where the uncertainty set is modeled with a polyhedral set) and describe the case study in Section 5.

Literature review

A review of major studies on the robust optimization and vaccine distribution SC models under pandemic conditions shows that reliable SCND facility locating models (with disorder confrontation approach) were first studied by [24] and, then, [25] conducted a review study on the robust SCND and categorized the uncertain nature of the SC parameters. During these years, there were other studies [26,27] that examined other SC features.

The current study has proposed a robust SCND model for a sustainable vaccine distribution that considers the economic, environmental and social dimensions. Here, the important point is how the model deals with uncertain parameters because their nature has further strengthened the use of robust data-driven optimization approaches discussed next.

Robust optimization

Robust optimization was first proposed by [28] in two dimensions: 1) optimality robustness which means the model will be close to optimal and 2) feasibility robustness which means the model will stay close to the feasible space in almost all cases of the realized solutions; "close" depends highly on how the modeler defines it. While stochastic approaches are impractical in uncertain environments where the distribution function of uncertain parameters is not known, robust optimization methods are quite practical because uncertain parameters can be modeled continuously or as discrete scenarios. Robust optimization models in the literature are classified as: 1) scenario-based, 2) fuzzy programming-based and 3) uncertainty interval-based.

In group 1, robustness measuring criteria minimize both the maximum cost and maximum regret (values farthest from the optimum). Some studies have minimized: 1) the maximum regret in the RSCND (closed-loop design) [29,30], 2) the average regret [31] and 3) the maximum cost [32,33], and some [34], [35], [36] have used constraints to ensure that the regret may not exceed a certain limit. Percent risk is another criterion used to measure the robustness of SC models [37], [38], [39].

In group 2, robust models are based on fuzzy sets. In SCNDs with uncertain parameters, the literature has repeatedly made use of the fuzzy programming based either on “flexibility” in operators, especially constraints, or “possibility” of the distribution of uncertain parameters. All studies (e.g., [41], [42], [43]) on these two methods have been solved in a general structure that converts fuzzy models to crisp models [40].

In group 3 that relates to uncertainty interval-based models, when uncertainty is associated with intervals, robust optimizations generally involve a high degree of conservatism; in other words, the problem infeasibility is highly protected. Study [44] was the first in this group that was aimed to present a model that kept the problem feasible for almost each feasible solution in a predefined set for all uncertain parameter realizations, but its problem was its high conservatism level. Studies [45] and [46] were proposed to reduce this level, but their approaches involved highly nonlinear formulations. For a better performance, Bertsimas [47] proposed a different approach where conservatism was controllable and, most importantly, its formulations were linear. While all previous models made real-world problem-decisions in only one step, subsequent studies [48] showed that they were usually made in several steps. “Distributionally robust optimization” is the title recent studies [49,50] have given to the uncertain distribution function in the robust optimization.

Pandemic vaccine supply chain

The vaccine SC in pandemic conditions has severe operational and time constraints. Uscher et al. [51] studied the influenza vaccine in 48 countries in 28 of which the high-risk vulnerable groups had the injection priority to reduce large-scale pandemic risks. To distribute the H1N1 vaccine, Araz et al. [52] divided the state of Arizona into 15 districts prioritizing 4 based on the estimates of the pandemic peak periods, population-based distribution sequence and a number of other factors.

Considering communications and ethics, Medlock and Galvani [53] developed an influenza vaccine distribution model that based allocations on such criteria as death, infection, age, and economic costs. Believing that the 6-12 year-old children were the cause of the inter-individual virus transmission because of their high activity/contact, in spring 2009 in Mexico, Lee et al. [54] considered the 20-39 age range eligible to receive H1N1 vaccine and developed a nonlinear dynamic model for its optimal allocation. Prioritizing the helpless and homeless for the H1N1 virus vaccination during a pandemic, Buccieri and Guetz [55] set up a case in Toronto, Canada, and vaccinated 38% of them in the existing clinics. Identifying the vaccine supply-distribution system as the most important factor, Davila-Payan et al. [56] developed a multivariate linear regression model to vaccinate the children and high-risk people against the H1N1 virus pandemic. Using a push-based strategy to distribute the vaccine in different areas, Huang at al. [57] conducted a study in Texas, USA, to categorize pregnant women, infants, and high-risk people as the priority groups. Chen et al. [58] developed an age-based simulation model to find the optimal COVID-19 vaccine allocation strategy by dividing the population into seven parts each with five age groups. They examined both dynamic and static policies and showed that first the older and then the younger groups had to be vaccinated.

An optimal vaccine SCND requires establishing an effective supplier-health department relationship [59]. In this regard, Marcello et al. [60] conducted a study in NY, USA, in 2009 on the vaccine distribution in an H1N1 pandemic using the data of medical centers to estimate the number of the distributed vaccines and showed that 96% of the suppliers were satisfied with the information recorded for the vaccine SC management. Covering about 89% of the drug stores in their study area for vaccine distribution, Fitzgerald et al. [61] conducted a study on the coordination/cooperation of suppliers, drug stores and public health centers across the US. In a study aimed to examine the vaccine supply-distribution methods used by governments in developed countries, Turner [62] showed that the traditional influenza vaccine distribution policies were weak and worrying tools and involved significant delays. To analyze different vaccination scenarios in the USA to prevent people from being hospitalized and die, Biggerstaff et al. [63] implemented a spreadsheet-base model the variables of which were the number of patients, different vaccination times, number of doses prescribed for various age groups and the hospitalization rate. They found that the starting time of the pandemic and vaccination were very influential factors in the vaccination program.

In a recent study, Govindan et al. [64] used the Fuzzy Inference System (FIS) for the demand management in healthcare systems to help reduce the devastating effects of COVID-19 disease by dividing the community into four groups based on the immunity system - age and underlying diseases - and showed that their classification helped fighting the disease in health systems. Tavana et al. [65] presented a mixed-integer linear programming model for equitable COVID-19 vaccine distribution in developing countries, Shim [66], Roy et al. [67] and Abbasi et al. [68] optimized the allocation of limited COVID-19 vaccine supply, Santini [69] assigned swabs and reagent for PCR testing, and other studies [70,71] considered the service level for herd immunity and capacity planning in the COVID-19 vaccine SC; however, none investigated the uncertain vaccine accessibility in developing countries. On the pandemic crisis, recently, several review studies have been conducted based on challenges of the COVID-19 pandemic situation in vaccine resilience supply chain and ripple effect in supply chains, including [72], [73], [74], and [75]. They categorized challenges in the literature and via vaccine supply chain (VSC) that the highest-importance are as follows:

Lack of vaccine monitoring bodies

Limited number of vaccine manufacturing companies

Unavailability of vaccines for developing countries

Lack of proper planning and scheduling

Vaccination cost and lack of financial support for vaccine purchase

Lack of correspondence between the VSC members

Considering supply chain challenges before a vaccine is administered to the general population can help vaccines successfully distribute vaccines. Therefore, solving critical challenges to the COVID-19 VSC is customary for a sustainable VSC [76] that could help the countries around the world to get out of the pandemic. So, research papers must follow to solve these challenges in the vaccine supply chain as this paper is organized to this goal.

Studies on the vaccine SC under worldwide pandemic conditions (supposed to continue till 2022), as well as the world vaccine suppliers, suggest that COVID-19 vaccine SCNDs are necessary, more than ever, to vaccinate global communities. But, supplying vaccines for less developed communities have an uncertain nature because the demand is high and suppliers are limited [77]. Here, robust models that yield robust solutions against uncertain vaccine-accessibility parameters seem appropriate because they do not need distribution functions for uncertain parameters; however, their high conservatism that imposes high costs on the system is a main weakness that may hinder their implementation, but being in the age of data and using the merits of robust data-driven optimization models can be a key to the solution of this problem.

As COVID-19 has long epidemic periods, sustainability is an issue worth considering in its vaccine SC, and since the literature is rather poor in robust studies that discuss the SC's economic, environmental and social dimensions, the present study has considered the mentioned dimensions and addressed the COVID-19 vaccine SCND at both the strategic level (construction site and warehouse capacity decisions) and operational level (inventory and product flow decisions). Its main innovations are as follows:

Addressing robust models’ high conservatism through a robust data-driven model.

Using cutting hyperplanes to implement more accurate and realistic data in the mentioned model.

Comparing the proposed DDRO approach with the box and polyhedral uncertainty set.

Supplying vaccine from several different domestic/foreign suppliers.

Addressing sustainability in three vaccine SC indicators by analyzing the environmental and social effects.

Implementing the social dimension in three indicators (No. of created jobs, days lost due to damage, potential hazards)

Implementing a case study in Iran to validate the proposed approach.

Problem statement

The proposed problem addresses the network design of an integrated producer-to-consumer SC that sends the vaccine to the end customer while packed and stored under its own specific conditions considering the vaccine flow and right supplier selection.

SCND of COVID-19 vaccine

One way to prevent infection to spread in the human body during a COVID-19 pandemic is vaccination by which communities can ensure long life, health and well-being (Pfizer, 2017), but since its outbreak is vast and multiple, the production, distribution and quality control of its vaccine have become a serious issue.

As the pandemic is now prevailing, minimizing its related complexities through a producer-distributor-consumer coordination can be effective in controlling the people's health [78]. The vaccine SC complexities are due to the production difficulties (manufacturing, mass production, etc.), globalization (global-level supply-demand balancing leading to adequate distribution), and observation of rules and regulations governing the development and manufacture of the vaccine. Since the international community needs public safety, the vaccine production and distribution trends will move towards higher quality and compliance with the related standards. Hence, to better respond to the complexities of the vaccine distribution network, it is necessary to optimize the planning and implementation of the SC network structure.

As production and distribution processes are two main vaccine SC features that make collective safety and health accessible to communities, they need special attention for design and planning. The SC in Fig. 1 includes domestic and foreign suppliers, international packaging centers, storage and distribution centers and provincial health centers as demand points. The chain, considering accessibility to a limited capacity of domestic vaccines, transfers them from the main producer in a foreign country to international packaging centers in the same country, then sends them to the destination country after packaging, keeps them in the storage and distribution centers in that country and, finally, transfers them to central health centers as points of demand by internal transportation.

Fig. 1

COVID-19 vaccine supply network.

COVID-19 vaccine can be stored in the production, packaging storage and distribution centers and the model determines its optimum storage level during the planning time horizon. The vaccine storage is possible by refrigerating and freezing and construction of storage and distribution centers can be based on either of them. The SC inter-level goods transportation is both internal (by road/railway) and external (sea/air). In the proposed SC, decisions are strategic and tactical (optimized by the model); the former determines the site/capacity of the storage/distribution/packaging centers and the related technologies and the latter include inventory, rate of product flow, supplier selection (in each time period) and the transportation mode. The proposed mathematical model has been so formulated as to minimize costs and environmental effects and maximize the positive effects of the SC activities as its social responsibility. It uses a robust optimization model to deal with the uncertain nature of supplying foreign vaccines for developing countries, but since its conservatism level is high, a DDRO-cutting planes model has been proposed for its solution. The measuring criteria of the social/environmental effects are examined in subsequent lines.

Environmental assessments

The COVID-19 vaccine SCND sustainability requires the effects of all its activities to be analyzed both environmentally and socially [79], [80], [81]; the former needs an accurate and standard measurement criterion from the highest to the lowest chain level. For this purpose, the life cycle assessment (LCA) is the most common structure founded based on ISO 14040 and 14044, but since its direct use is both costly and time consuming, use is made of the “ReCiPe” method that calculates the environmental effects based on such activities as transportation, producing various materials and inventory maintenance using the SimaPro commercial software and an Eco-indicator 99 Standard-based updating database [82]; the same have been used in the current study to assess the environmental effects.

Social assessments

Measuring the social responsibility is generally difficult because it has a complex nature and lacks a specific criterion and depends on several stakeholders and disciplines each with its own size and criterion. Eventually, as organizations and companies considered the sustainable SC development an important issue, the International Standard Organization (ISO) recently introduced the ISO2600 measurement standard which is a measuring criterion for the social responsibility and is classified into: 1) organizational governance, 2) human rights, 3) labor practices, 4) environment, 5) fair operating practices, 6) consumer issues and 7) community involvement and development.

This study has tried to accurately estimate the chain's social impacts by identifying the related criteria in the COVID-19 vaccine production/distribution SC and use them in the mathematical modeling. According to ISO2600, stakeholder selection is of special importance as it is the first step in identifying these criteria. Table 1 that lists the concerns for categorizing these criteria has three social responsibility measurement criteria to find the social effects of the COVID-19 vaccine SC, which are integrated by weighting coefficients that not only point to the importance of the criteria, but are also used most to help integrate different criteria with different measurement units.

Table 1

Criteria for social responsibility.

ISO 26000 categories	evaluation criteria	Social effects	Shareholder
Labor practices	Days lost due to damages	Job security	workers
	Number of created job opportunities	Occupational justice
The environment	Contaminating products	Environmental vulnerability	Society

Following the Covid-19 pandemic and losses of different jobs worldwide, developing novel SCs for its vaccines can effectively create new jobs to avoid its adverse social effects, but the pandemic's temporary nature makes the job security a critical issue for the employees. Expert human resources working in vaccine production plants are usually exposed to being contaminated with chemical products; these are issues necessary to consider for the SCs’ social sustainability.

Unique features of the proposed vaccine SC model that make it different from other models are: 1) establishing separate distribution and packaging centers for special vaccine maintenance, 2) showing the importance of the social dimension of a sustainable development, especially in the COVID-19 pandemic situation, 3) paying attention to different transportation modes in a suitable substrate for all routes for special vaccine maintenance, 4) considering uncertainty of the foreign vaccine accessibility (because of inter-buyers competitions), 5) paying attention to supplier selection based on several COVID vaccine suppliers all over the world and 6) supplying vaccines from domestic and foreign suppliers.

Mathematical modeling

Notations of the proposed mathematical model are provided in the Appendix and its objective functions are introduced next to develop its structural constraints.

Objective functions

Based on the sustainability concept, the objective functions considered for COVID-19 vaccine SCND are economic, environmental and social.where is the sum of the fixed costs of constructing: 1) international packaging centers for delivery to other countries, 2) distribution centers in the destination country and 3) domestic production centers.where is the variable costs of providing capacity in each SC facility; internationally, the capacity involves those of the packaging centers, domestic distribution centers and domestic producers in each period.where is the inventory costs for each facility for, respectively, those in the international manufacturers and packaging centers, and those in the domestic producers and distribution centers in the destination country in each period of the planning horizon. where , and Eqs. (4), (5) and (6) show, respectively, the vaccine production, distribution and packaging costs in the related international and domestic centers in each period of the planning horizon.where shows the shortage-related costs in each international production, packaging and distribution center and in domestic production centers in each period of the planning horizon.where is the capacity expansion cost assignable to each facility if the design requires the expansion of international production, packaging and distribution centers and domestic production centers, respectively.where (a prominent part of the economic objective function of costs associated with transportation) relates, respectively, to the costs of transportation from the international supplier to packaging centers, then to distribution centers in the destination country, then to provincial medical centers and also from domestic suppliers to medical centers in each period of the planning horizon. where (Eq. 10) is the supply costs of raw materials needed to produce the vaccine in both international and domestic centers in each time period and (Eq. 11) is the fixed costs of purchasing road vehicles.where is the sum of the periodic costs.which is the final COVID-19 vaccine SC economic objective function where Z is minimized with the interest rate (IR) parameter and the net present value relation plus time-independent costs minus road vehicles’ scrap costs.where is the environmental emissions due to vaccine transportation from the foreign supplier to packaging centers, then to distribution centers, then to medical centers and, finally, from the domestic producer to medical centers in each time period. where , EFFR and (Eqs. 15, 16, 17, respectively) are the environmental effects due to the production process in both domestic and foreign plants (Eq. 15), the construction of international packaging facilities, distribution centers and domestic production centers (Eq. 16) and the packaging process in international centers (Eq. 17) in each time period.where is the environmental effects of the vaccine maintenance in international production, packaging and distribution centers and domestic production centers in each time period.where is the sum of the environmental effects of the time-dependent variables.which is the final environmental objective function.where is the number of jobs created by transportation from the foreign manufacturer to the packaging center, then to the domestic storage and distribution centers, and then to provincial health centers, and also among domestic manufacturers and health centers.where is the number of jobs created in each domestic/foreign vaccine production center. where (Eq. 23) is the number of jobs created based on the activities of international packaging centers and (Eq. 24) is that created in domestic distribution centers.where is the first term of the social objective function and shows the sum of the jobs created in all SC processes.where is the second term of the social objective function and is related to the people's number of lost days due to injuries in the workplace.where is the last term of the social objective function and is related to the workplace risks calculated only in the production process in foreign and domestic centers.

Components of the economic objective function:

Components of the environmental objective function:

Components of the social objective function:

Now, all the current SC social effects are formulated as follows by giving weight to each component; maximizing number of jobs due to SC activities with weight (with +ive sign), number of days lost due to injuries with weight (with -ive sign for minimization) and environmental risks with a weight (with -ive sign for minimization):

Constraints

This section presents the model constraints on the capacity (facility and inventory), transportation, construction, and capacity expansion and set of constraints on vaccine flow and inventory balance.

Capacity constraints

This section deals with the constraints on the capacity to ensure that the capacity of each domestic/foreign production center is less than or equal to its designed capacity in each period; this also applies to international packaging and domestic distribution centers that depend on different storage technologies (29), (30), (31), (32) formulate these constraints for, respectively, foreign production centers, domestic production centers and packaging centers.

The above formulae limit the production of each facility based on its designed capacity, but here, effort is made to rationally determine this capacity by determining the related upper and lower limits in the construction and economic dimensions. Eq. (33) specifies the design capacity of a foreign production center () with two upper and lower limits (); this trend is repeated for domestic production centers, foreign packaging centers and distribution centers in (34), (35), (36).

Here, capacity-related constraints are considered for centers capable of keeping inventory. Eq. (37) limits the foreign supplier inventory () by maximizing the amount storable in that center (). These constraints have also been formulated for domestic packaging, distribution and production centers in Eqs. (38)-(40).

Constraints related to the flow of raw materials to domestic and foreign production centers are next. It should be ensured that the raw material required for vaccine production should be less than or equal to the maximum amount available in each time period. Eq. (41) is related to the raw material supply for foreign producers and Eq. (42) is related to that for domestic producers.

Capacity expansion constraints

Capacity expansion is an important concept aimed to create a progressive capability for the capacity of the designed facility to be available when needed to better meet the demand. To formulate this concept, use is made of the designed capacity in each period which equals that of the previous period plus the capacity expansion in that period. This definition has been formulated for all four types of facilities (domestic and foreign production centers, packaging centers and distribution centers) in Eqs. (43)-(46).

Transportation constraints

As mentioned before, international transportation is by sea and air and the domestic transportation is by railway and road. Constraints related the capacity of the transportation fleet are () and those on accessibility to a sufficient number of vehicles are () Eqs. (47) and (48) relate to the vaccine transportation within each country from production centers to packaging centers and from distribution centers to provincial medical centers, Eq. (49) is used for the same fleet for domestically produced vaccine delivery and Eq. (50) is for international transportation by sea and air.

Flow balance constraints

Constraints discussed here try to balance the vaccine flow in different COVID-19 vaccine SC layers; formulas (51-55) try to match the demand and supply by balancing the flow and correct SC modeling. Eq. (51) uses the -type rate () to calculate the vaccine production in foreign centers based on the available raw materials, Eq. (52) uses the vaccine packaging coefficient () to determine the packages prepared for external shipment, Eq. (53) uses the distribution coefficient to convert part of the shipment to distribution packages () after the international shipment and arrival of packages in the destination country and Eq. (54) uses the relevant factor () to prepare the raw materials for vaccine production in a conversion process in domestic centers. Obviously, vaccines produced and sent from both foreign and domestic production centers to provincial medical centers should meet the demand of each center in each period (), but since this parameter is uncertain, the accessibility of each vaccine in the competitive world market (), is calculated in Eq. (55).

Inventory constraints

Eqs. (56)-(59) have been formulated to create balance in centers that keep inventory of goods. These values are calculated assuming that foreign and domestic production, packaging and distribution centers store inventories. Eq. (56) shows the inventory balance in foreign production centers, which should be equal to that center's outflow, that is, production in that period plus the reduced inventory minus the shortage Eqs. (57)-(59) show the same concept for the packaging, distribution and domestic production centers.

Other constraints

Eq. (60) shows the rational constraint on the selection of a vaccine storage technology in a distribution center and Eq. (61) re-determines the status of each variable.

Data-driven robust optimization (DDRO)

In dealing with uncertain parameters, the robust optimization has been generally popular with researchers in recent years and used more than previous uncertain stochastic programming approaches. This tendency has been greater in SCNDs due to high investments and information inaccessibility at the beginning of the design. Selecting robust optimization approaches is due, mainly, to their independence from the distribution function of uncertain parameters because they can provide an accurate estimate of such parameters by only knowing their variation interval; access to their distribution is not necessary.

The robust optimization method proposed in this study has higher conservatism than its rivals; it defines a closed convex set for the problem's uncertain parameters and presents, on its basis, a closed formula of its robust counterpart. Here, the COVID-19 vaccine SC is executed by first introducing and formulating a robust optimization method based on closed uncertainty sets and then presenting the proposed approach.

Robust convex optimization

To see how to address set-oriented robust uncertain parameters, consider the following linear problem: where and are uncertain and their variations, in symmetrical and entirely random intervals, are and , respectively; and show their nominal values and and are their variation rates in each interval. An uncertain parameter can be rewritten as follows:where is a random variable belonging to set that has a symmetric distribution in the range [-1,1]. If this variable is multiplied by the variation rate of each uncertain parameter interval and its nominal value is added, the result will be a mathematical representation of the uncertain parameter based on its changes. Substituting the closed form obtained above in Eq. (63) for and will yield:where and belong to polyhedral set (Fig. 2 ) where for the values on the right. Using a random variable in the form “”, this set is defined with control parameter as follows to limit the uncertainty space:

Fig. 2

Polyhedral uncertainty set.

Considering the definitions, Boyd et al. (2004) [83] presented a linear robust twin with a high degree of conservatism as follows for the maximization problem in constraint (65):

Applying cutting hyperplanes to control conservatism

The approach discussed in the previous section belongs to the category of static robust optimization methods that have a fixed routine throughout the planning horizon. The main weakness of such methods is their high conservatism and paying special attention to the occurrence of the worst case of the uncertain parameter. In SCNDs, this weakness imposes high costs on the system mostly because of lying in a safe point preventing the occurrence of the predicted state with high probability.

To address the high conservatism of static robust optimization approaches, researchers have, in recent years, developed new dynamic robust optimization methods that highly reduce conservatism in classical models by realizing some data of the uncertain parameters over time and limiting, on this basis, their set [84]. Fig. 3 shows the decision-making structure in robust dynamic programming by a numerical example.

Fig. 3

structure of decision making in dynamic robust programming.

Dynamism in the DDRO modeling based on closed convex sets reduces the conservatism level. Among conservatism reduction approaches, optimization methods, e.g. cutting hyperplanes, limit the space of the set of uncertain parameters based on the data realization [85]. After the realization of some of the data, these planes are plotted on the set of uncertain parameters, in the next time period, based on these data (Fig. 4 ) to not only limit its high conservatism, but also improve the solution robustness.

Fig. 4

Polyhedral uncertainty set with cutting hyperplanes.

Given the number of cutting hyperplanes (C) with gradient vector and intercept, hyperplane is defined with random vector as follows: where is the number of uncertain parameters in the i h constraint.

The considered cutting planes consist of a coefficient () and a constant () the equation of which is formulated by a random variable. In this study, the DDRO approach has been implemented based on a polyhedral packaging uncertainty set where cutting planes reduce conservatism and increase robustness by limiting these areas. Next, effort is made to obtain a linear closed form of the common area between the initial sets and the cutting planes.the inner maximization equation of which is as follows: where and .

Hence, the dual variable is and the dual problem is as follows: simplified as:

As, the closed form of Eq. (71) is as follows:

The final linear form of which with auxiliary variables and is:

Formation of cutting planes

To keep the convexity of the uncertainty space, use is made of linear cutting planes with the following general form:where is the coefficients’ linear matrix and is a vector of constant values; accordingly, the linear cutting planes will be where is the number of cutting planes and is the cardinality of the set of uncertain parameters. The optimization model of achieving the components of these planes is as follows: where is a sample of uncertainty space , specifies its size, is the sample distance to cutting plane, , is a normal vector with states to create different cutting planes and is a threshold to be specified each time. The final equation of cutting planes is as follows:

Case study

In Dec. 2019, authorities in Wuhan City, Hubei Province, China, reported some unknown virus-caused pneumonia cases related to a seafood wholesale market [86] which was first referred to as a Mysterious Respiratory Disease in China and then renamed as COVID-19 by the WHO [87]. COVID-19 is caused by an emerging bat-origin corona virus (SARS-CoV-2) not previously reported in humans [88]. It is a highly acute and contagious disease began in China in mid-December and spread rapidly during the New Year holidays because information was none and trips were not restricted.

Prioritizing prevention over treatment has long been an issue in the world health system and the sudden onset of the COVID-19 crisis reaffirmed this priority. Inability to control this disease has put much pressure on the countries’ treatment systems and consequences of the related imposed costs have had such profound effects on their economies that early vaccine accessibility has become their first priority to appropriately respond to and prevent the crisis. Vaccines, approved by the WHO till April 2021, include Sputnik V, Pfizer, Astrazeneka, Modrena and Sinopharm. COVIran Barakat is an Iranian vaccine that has passed its third clinical phase, entered the treatment system for general vaccination and used first to vaccinate the leader of Iran to show his respect for young Iranian scientists.

The COVID-19 vaccine SC consists of such complicated systems as processing, equipment, means of transportation and afflicted places that need proper temperature for their inventory. This SC is aimed mainly to provide all humans with a better vaccine access and distribution and since failure to understand its components and interactions can make the vaccine accessibility difficult around the world, a robust SCND can further enhance the vaccine accessibility, especially in developing countries. The present study has implemented its model in a practical case in Iran.

With 2,823,887 affected and 77,994 deaths as of May 2021 [89], Iran has been one of the twenty countries with the highest COVID-19 afflicted and present in all the related peaks. Various disease outbreaks made the government try to not only control it with such policies as closing various centers nationwide and enacting quarantine-related laws at different intervals, but also consider a comprehensive nationwide vaccination plan on its agenda. The present case study has addressed the Russian Sputnik V, Indian Astrazeneka, and Chinese Sinopharm vaccines ordered by the government of Iran, and also the Iranian COVIran Barakat that has passed its third clinical phase. It has entered the treatment system to supply part of the vaccine needed for a general nationwide vaccination (there are still other domestic vaccines underway to be dealt with later).

Data collection for the COVID-19 vaccine SCND

This section presents the data-collection sources for the input parameters of the COVID-19 vaccine SC to fruitfully conduct the case study by obtaining valid information. As the temperature for the vaccine storage/transportation was — 80° C, the fleet used in this study operated accordingly. The information related to the foreign transportation, including costs and accesses in Russia, China and India, was gathered from the Ministry of Roads and Urban Development of Iran [90], that of the international shipping was collected from Iran Chamber of Commerce [91], that on the size of the vaccine packages shipped to Iran was provided from the documents passed among the Ministry of Health of Iran and those of the three mentioned countries and that on air/sea transportation was provided from the Civil Aviation Organization of Iran [92] and the Shipping Organization of Iran [93], respectively Table 2 [94] lists the details, characteristics and information of vaccines used in this study including the platform, number of doses, schedule of doses, route of injection and its developer.

Table 2

Information on international COVID-19 vaccines.

Developers	Route of administration	Schedule	Number of doses	Platform	Vaccine
"Gamaleya Research Institute; Health Ministry of the Russian Federation"	IM	Day 0 + 21	2	Viral vector (Non-replicating)	Sputnik V
Sinopharm + China National Biotec Group Co + Wuhan Institute of Biological Products	IM	Day 0 + 21	2	Inactivated virus	Sinopharm
AstraZeneca + University of Oxford	IM	Day 0 + 28	2	Viral vector (Non-replicating)	AstraZeneca
Barakat Foundation (Execution of Imam Khomeini's Order)	IM	Day 0 + 28	2	Inactivated virus	COVIran Barakat

As the mathematical model was to locate such facilities as the foreign packaging and domestic distribution and production centers, potential sites for the latter were considered at two provinces: 1) Razi Vaccine and Serum Research Institute (Alborz Province) and 2) Pasteur Institute of Iran (Tehran Province) considering accessibility and available resources/facilities. Candidate provinces, in Iran, for distribution centers were Tehran, Tabriz, Markazi, Fars, Isfahan, Khuzestan, Khorasan Razavi, Kermanshah and Hormozgan and potential packing centers were Beijing, Shanghai and Guangzhou in China, Moscow, St. Petersburg and Kazan in Russia and Mumbai and Delhi in India; all provinces in Iran were considered for provincial medical centers.

This final section is aimed to address the vaccine demand and injection priorities. According to the strategic vaccination document approved by the Iranian government, groups prioritized for vaccination are (Fig. 5 ): 1) medical staff, 2) the elderly and those with underlying diseases, 3) those with sensitive jobs and high public contacts and 4) the public; the related information has been collected from the “Ministry of Health” and “Statistics Center” of Iran [89,95]. It is worth noting that since all vaccines are injected in two doses, demands in periods after the first will be the total of the first dose for new people and the second dose for those already received the first injection. The model presented in this study considers a 12-month planning time horizon for this project.

Fig. 5

Vaccination phasing.

In consultation with Iranian social pathology experts, weight coefficients were determined for critical social effects in the Covid-19 vaccine SC resulting in 0.41 for occupational justice, 0.3 for environmental vulnerability and, finally, 0.29 for job security.

Implementation of the COVID-19 vaccine SC

The proposed multi-objective linear mathematical optimization model with three conflicting objective functions that face uncertainty in the vaccine accessibility parameter, has been solved with the GAMS 24.1.2 software using a quad-core, 2GHz frequency, 8 GB RAM processor and the epsilon constraint approach. Among the set of obtained Pareto solutions, the authors have selected the one that yields efficient and same-level values ​​in all objectives.

The linear cutting planes, consisting of coefficients and fixed values, were formed by the mathematical DDRO model implemented based on closed uncertainty sets using data sampling and models presented in Eqs. (76)-(77). As the foreign vaccine accessibility parameter () (3 in this study) is uncertain, planes are required (Table 3 ) to proceed. Applying these planes on the multidimensional uncertain parameter set yielded the final model results.

Table 3

Cutting hyperplanes.

Cuttingplane:h_c(ξ)=q_cjξ_ij+d_c
0.3666ξ_2+0.7074=0	c=5	0.4521ξ_1+0.6666=0	c=1
0.4258ξ_1+0.1295ξ_3+0.8148=0	c=6	0.5332ξ_1+0.3604ξ_2+0.8937=0	c=2
0.2489ξ_3+0.7702=0	c=7	0.2322ξ_1+0.65=0	c=3
0.4772ξ_1+0.6233ξ_2+0.8747ξ_3+0.6141=0	c=8	0.2318ξ_2+0.4058ξ_3+0.8694=0	c=4

Next, the SC model was implemented after applying the cutting planes formed in the set of uncertain accessibility parameters of each vaccine type Fig. 6 shows the SC network that supplied the COVID-19 vaccine for Iran in an international dimension. Packaging centers for distribution purposes were selected, from among the potential locations, to be in Moscow, Beijing and Delhi based on the mathematical model solution results (), and 8000 of the Russian, 9000 of the Chinese and 10000 of the Indian vaccine doses were considered in each package for international shipments.

Fig. 6

International supply network.

After the international packaging centers were determined and the optimal vaccine flow outside Iran found its final form, the SC form was addressed inside the country. Out of nine potential locations for foreign-to-domestic vaccine distribution centers, the optimal solution selected only 5 - Tehran, Isfahan, Ahvaz, Kermanshah and Bandar Abbas - to receive the vaccines bought from foreign manufacturers and ship them to the desired distribution center (); on the other side of the defined SC lay the domestic (COVIran Barakat vaccine) production centers: 1) Pasteur Institute of Iran, in Tehran and 2) Razi Vaccine and Serum Research Institute, in Alborz (Fig. 7 ) (). Regardless of the type, all four available vaccines types - three imported and one domestic - can be used to meet the nationwide demand.

Fig. 7

Domestic SC network.

The previous section located the construction sites of 5 distribution centers of the imported vaccines all of which used all three vaccine types to meet the demand of the study population. Table 4 shows the rate of the vaccine distribution flow to provincial medical centers in each period for different centers and vaccine type (). As shown, since the vaccine production and transportation costs are different, the model has used some vaccines more in certain periods to reduce costs. Table 5 lists the vaccines produced in domestic centers (Fig. 7) in each time period separately for the related production center (). Tables (4) and (5) show an appropriate foreign-domestic production balance in the optimal solution. It is obvious in Table (5) that the domestic vaccine accessibility in the early period is far from expectation considering the commencement of the production process. Regarding road transportation, the model has determined the need for 108 truck units for all SC displacements.

Table 4

International vaccine flow to provincial medical centers.

		Time Horizon
Distributions centers	Vaccine Type	t=1	t=2	t=3	t=4	t=5	t=6	t=7	t=8	t=9	t=10	t=11	t=12
Tehran	Spontik V	1200927	1355390	1554804	1554804	1045778	-	1157197	1554804	1554804	1200927	1014488	1554804
	Sinopharm	321510	305100	309800	309800	-	-	309800	309800	309800	309800	390056	824126
	AstraZeneca	515900	515900	291454	515900	-	-	636800	636800	636800	636800	636800	1788600
Esfahan	Spontik V	301755	426000	426000	393394	-	-	234707	239561	426000	258389	-	-
	Sinopharm	82890	268110	351000	84800	80100	351000	351000	80100	351000	351000	270743	439074
	AstraZeneca	1964470	1622155	2186800	2186800	1539975	-	2186800	2186800	2186800	2186800	2186800	2186800
Ahvaz	Spontik V	446000	446000	-	-	-	-	-	-	-	-	-	-
	Sinopharm	1189926	1339200	1339200	1339200	-	-	1339200	1339200	1339200	1339200	1339200	1339200
	AstraZeneca	327195	445500	445500	2181300	-	-	2181300	2181300	2181300	2181300	2181300	2181300
Kermanshah	Spontik V	218589	312325	411000	411000	-	-	61976	296683	411000	411000	-	-
	Sinopharm	22535	237330	270900	270900	-	-	270900	270900	-	-	-	-
	AstraZeneca	19931	928429	1920600	1920600	-	-	1920600	1920600	1920600	1920600	1920600	1920600
Bandar-abbas	Spontik V	1463661	1986000	2283900	1986000	1986000	1986000	1986000	1986000	1986000	1986000	1986000	1986000
	Sinopharm	326700	-	-	-	-	-	326700	-	-	-	-	-
	AstraZeneca	1303999	189726	1706624	1974500	1974500	1974500	1974500	1974500	1974500	1974500	1974500	1974500

Table 5

Domestic vaccine flow to provincial medical centers.

	Time Horizon
Manufacturer	t=1	t=2	t=3	t=4	t=5	t=6	t=7	t=8	t=9	t=10	t=11	t=12
Razi technology Incubator	-	2249429	804219	2435860	2238037	-	-	2943422	2943422	2943422	2943422	4143422
Pasteur institute of Iran	-	3793579	4658872	2956487	1665899	613472	281078	4671100	3000000	3000000	3000000	6815479

Fig. 8 shows which foreign or domestic vaccine producer covers the provincial medical centers to meet the demand for a nationwide vaccination. Supplying by internal vaccines shows that not all medical centers have been supplied with even a minimum amount due, maybe, to the limited domestic production capacity in the study period. In the view related to foreign vaccines, each medical center is often covered by several distributors reducing the transportation costs.

Fig. 8

Demand Centers’ Allocation Network.

Meeting the demand for COVID-19 vaccine through both domestic and foreign producers requires careful planning of their capacities. Fig. 9 shows the production capacities () of the Russian Sputnik V, Chinese Sinopharm and Indian Astrazeneka designed for uncertain parameters of accessibility and ancillary costs with the maximum capacity relating first to the Russian and then to the Indian vaccines to meet the demand in the early years. Fig. 10 shows the design capacity variable () for domestic production centers where the capacity specified for “Pasteur” is somewhat more than that for “Razi”. This can also be concluded from Fig. 8 where the coverage of this center is more than that of the Razi Institute.

Fig. 9

Capacity of international vaccine manufacturers.

Fig. 10

Capacity of domestic vaccine manufacturers.

Next, the proposed mathematical DDRO model is compared with the deterministic programming and box and polyhedral uncertainty set-based approaches to evaluate its performance.

DDRO model performance

The rapid corona disease epidemic boosted the demands for its vaccine in developed countries after it was produced and made its access uncertain for developing countries. The proposed data-driven approach has addressed this uncertainty and is validated next for each foreign vaccine's accessibility.

To this end results of the proposed DDRO approach are considered as the solutions of the mathematical model and compared with the box and polyhedral set in two uncertainty modes for the accessibility parameter based on the classical convex programming; in other words, solutions of other methods are compared with those of the proposed approach to explain its merits. To observe the differences between different approaches and show the proposed model performance better, the uncertainty interval of the respective parameter,, has been changed ( is the distance between the real and realized data).

Fig. (11) shows the comparisons made in the present study based on the economic cost objective function. Performance of the three approaches facing uncertainty is shown by changing the real-realized data difference parameter in its uncertainty interval. In the box uncertainty set-based approach, solutions are much worse than others due to the method's high conservatism that ensures solutions have high robustness against the variations of the uncertain parameter, but in the polyhedral uncertainty set-based approach, the objective function shows lower costs and solutions are improved and have less conservatism, but still impose high costs on the system for cases that are less likely to occur. In short, the proposed DDRO method greatly reduces the conservatism of classical, uncertainty set-based approaches and the solutions, while quite robust, are much better than those of other approaches.

Fig. (12) compares the solutions in different approaches based on the economic objective function and the difference between the real and realized values. As shown, the classical box/polyhedral uncertainty set-based approaches have created high conservatism in solutions, causing the environmental objective function to get worse, but the proposed approach has shown acceptable performance in the environmental dimension, reduced the high conservatism of classical approaches significantly and yielded a much more robust solution. This event is clearly obvious in Fig (13) in the behavior of the three-part social objective function.

Discussion and main results implications

In the vaccines distribution context, a pandemic VSC is different than that of a traditional VSC because governments are directly procuring vaccines from the manufacturers bypassing the traditional chains of wholesalers and distributors. Hence, healthcare experts and VSC analysts are looking for proper policies and adequate strategies for appropriate vaccine manufacturing and distribution to fight against the COVID-19 pandemic. Investigating the existing literature shows the vaccine supply chain has not been updated according to the pandemic situation. The main challenges that are not considered in COVID-19 vaccine supply chain literature and can affect supplying vaccines performance are as follow:

Integrating import and manufacture decisions simultaneously in countries with vaccine production knowledge

Planning vaccine supply chain in a sustainable form (especially the social aspect)

Addressing Covid-19 vaccine access uncertainty in the competitive world market

In this regard, the present study has been aimed to fill these gaps so that governments can remove the barriers to the public access to appropriate vaccines and according to the results there is considerably good integration to solve the problems. As the proposed study has been implemented in a real case study in Iran, so far among the first ten most affected countries, the managerial insights have been discussed for the results of the study in each section.

Simultaneous import-production decision integration

For countries that have the vaccine production knowledge, the proposed model enables an integrated approach that illustrates, in its network, the levels associated with domestic and foreign suppliers for the vaccine production and importation, respectively. In the COVID-19 vaccine issue, governments, which are generally the upstream-level decision makers, can benefit from the integrated decision-making by precise vaccine production plans. Benefits of this concurrence are: 1) reduced costs of supplying the Covid-19 vaccine, especially in a situation where the world is facing an economic crisis caused by the pandemic, 2) careful planning of the vaccine importation based on the production capacity of each case due to the high complexity of the vaccine importation supply chain, 3) establishing proper coordination with the world's main producers of the Covid-19 vaccine - Sputnik, Sinopharm and AstraZeneca - and 4) precise production planning for countries that have the vaccine production knowledge.

Using a real case study in Iran, the proposed model depicts the relations among the main foreign vaccine producers for importation in Fig. 6, locating the domestic production facilities and foreign vaccine distribution centers in Fig. 7, and allocating vaccine producers to meet the demand in provincial medical centers (Fig. 8) Tables 4 and 5 specify the optimal flow rate of the domestic and foreign vaccines based on the production and importation capacities of each vaccine over the horizon, and formulate the economic objective function to minimize costs and optimally balance the importation and production.

Sustainable planning of the vaccine supply chain

The environmental and social activists are nowadays seriously concerned about how wrongly the humans behave with the nature and believe it is a threat for the next generation's access to a healthy environment and community. This situation has worsened with the COVID-19 disease pandemic, loss of many jobs and increased environmental effects.

The proposed COVID-19 vaccine model considers the environmental and social dimensions of a sustainable SC network and formulates its second and third objective functions accordingly. The second objective function minimizes environmental effects of the entire network processes and ensures the related activists that the emission of greenhouse gases is controlled to a desirable level based on the ReCipe Standard. The third objective function addresses such social concerns as the number of created job opportunities, days lost due to the damage, and contaminating products in a situation where the manpower has lost many jobs due to the pandemic throughout the world. A comprehensive consideration of this dimension of the sustainable SC network based on the ISO2600 Standard ensures the human resources activists that, under the pandemic, the desired effects of the Job security, occupational justice and environmental vulnerability on the network's manpower will be maximized; results of these objective functions are shown in Figs. 12 and 13.

Fig. 12

Data-driven model performance versus other approaches on environmental function.

Fig. 13

Data-driven model performance versus other approaches on social function.

Uncertain COVID-19 vaccine accessibility in the competitive world market

A challenge facing a vaccine SC in a disease pandemic is the competitive markets where governments around the world try to provide the vaccine. Such markets can be due to the governments’ varying behavior such as providing beyond their needs, production limitations in producing countries, and so on that lead to uncertainties to access and purchase the vaccine.

This study has used modern, robust, data-driven optimization methods and implemented the model under uncertain conditions using cutting planes to present a realistic view of the uncertain nature of the vaccine accessibility in pandemic situations because such planes are extracted by realizing the vaccine accessibility data for purchasing. By controlling the conservatism, the proposed approach allows Govt. managers to make vaccine production-/importation-related decisions that conform more to reality and avoid high-conservatism decisions that impose high costs and adverse environmental/social effects on the natural system. Results obtained from this merit of the proposed approach are shown in Figs. 11, 12 and 13 in each of the related objective functions.

Fig. 11

Data-driven model performance versus other approaches on cost function.

Conclusions

The present study has addressed a robust COVID-19 vaccine SCND using a mathematical multi-period mix integer linear programming model implemented for a 1-year time horizon for a nationwide vaccination in Iran. It is a robust SC network with economic, environmental and social dimensions; the latter considers the jobs created and job days lost by workers due to risks and injuries in the workplace. Against the worldwide vaccine-access injustice, a robust data-driven model has been developed to address the uncertainty of the vaccine accessibility based on closed convex uncertainty sets.

A merit of the proposed approach is a 21% reduction in the conservatism of the classical robust box/polyhedral set-based approaches while solutions are robust and resilient against variations of uncertain parameters. In the strategic dimension too, the results indicate the construction of two domestic vaccine production centers in Razi and Pasteur institutes and five international vaccine distribution centers in Iran - Tehran, Isfahan, Ahvaz, Kermanshah and Bandar Abbas. The proposed model determined what foreign or domestic vaccine producers covered the provincial medical centers to meet the demand for a nationwide vaccination and optimized their designed capacities under uncertain vaccine accessibility based on the Iranian government's vaccine injection plan. This study can help governments to plan for resources, facilities, and distribution vaccines under competitive purchase situations.

Considering the further development of different types of vaccines in the world, future studies can: 1) focus on wider vaccine SC networks, 2) Re-modify the SC structure of such vaccines as Moderna and Pfizer considering different maintenance conditions and 3) Propose a new objective function to minimize the lead time to satisfy the demand more effectively.

Author Contribution Statement

The authors confirm contribution to the paper as follows:

Study conception and design: H. Gilani, H. Sahebi; data collection: H. Gilani; analysis and interpretation of results: H. Gilani, H. Sahebi; draft manuscript preparation: H. Gilani, final revision: H. Sahebi.

All authors reviewed the results and approved the final version of the manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix.

Index
I	Set of locations foreign vaccine manufactures
J	Set of potential locations for packaging centers
K	Set of potential locations distribution centers
R	Set of locations for provincial medical centers
V	Set of potential locations internal vaccine manufactures
L	Set of transportation mods
P	Set of vaccine types
M	Set of vaccine inventory technologies
T	Set of time periods
Technical parameters
FIXP_j	Fixed cost of packaging center at location i
FIXD_km	Fixed cost of distribution center at location k with inventory technology m
FIXI_v	Fixed cost of internal vaccine manufacturing center at location v
VCMipt	Variable cost of foreign manufacturing center at location i with vaccine type p in period t
VCPjpt	Variable cost of packaging center at location j with vaccine type p in period t
VCDkpmt	Variable cost of distribution center at location k with vaccine type p and inventory tech m in period t
VCIvt	Variable cost of internal vaccine manufacturing center at location v in period t
PCMipt	Unit production cost of vaccine type p at foreign manufacturing center i in period t
PCIvt	Unit production cost of internal vaccine at internal manufacturing center v in period t
PCDkpmt	Unit pre-distribution cost of foreign vaccine type p at distribution center k with inventory tech m in period t
PCDIvt	Unit pre-distribution cost of internal vaccine at internal manufacture center v in period t
PCPjpt	Unit packaging cost foreign vaccine type p at packaging centers j in period t
ICMipt	Inventory holding cost per unit of foreign vaccine type p at foreign vaccine manufacturing center i in period t
ICPjpt	Inventory holding cost per unit of foreign vaccine type p at packaging center j in period t
ICDkpmt	Inventory holding cost per unit of foreign vaccine type p at distribution center k with inventory tech m in period t
ICIvt	Inventory holding cost per unit of internal vaccine at internal vaccine manufacturing center v in period t
SCMipt	shortage cost per unit of foreign vaccine type p at foreign vaccine manufacturing center i in period t
SCPjpt	shortage cost per unit of foreign vaccine type p at packaging center j in period t
SCDkpmt	Shortage cost per unit of foreign vaccine type p at distribution center k with inventory tech m in period t
SCIvt	Shortage cost per unit of internal vaccine at internal vaccine manufacturing center v in period t
CTRVijpl	Unit transportation cost of foreign vaccine type p form foreign vaccine manufacturing center i to packaging center j by transportation mod l
CTRPjkpl	Unit transportation cost of foreign vaccine type p form packaging center j to distribution center k by transportation mod l
CTRHkrpl	Unit transportation cost of foreign vaccine type p form distribution center k to provincial medical center r by transportation mod l
CTRIvrl	Unit transportation cost of internal vaccine form internal vaccine manufacturing center v to provincial medical center r by transportation mod l
CCEMipt	Unit cost of capacity expansion of foreign vaccine type p at foreign vaccine manufacturing center i in period t
CCEPjpt	Unit cost of capacity expansion of foreign vaccine type p at packaging center j in period t
CCEDkpmt	Unit cost of capacity expansion of foreign vaccine type p at distribution center k with inventory tech m in period t
CCEIvt	Unit cost of capacity expansion of internal vaccine at internal vaccine manufacturing center v in period t
MV_l	Market value of vehicle with mod l (just trucks)
IR	Interest rate
nt	Number of periods
SUMAXpt	Maximum available raw material can be ordered for foreign vaccine type p in period t
SUIMX^t	Maximum available raw material can be ordered for internal vaccine in period t
PRNipt	Unit purchasing cost of raw material foreign vaccine type p by foreign vaccine manufacturing center i in period t
PRIvt	Unit purchasing cost of raw material internal vaccine by internal vaccine manufacturing center v in period t
INR_l	Unit fix purchase cost of transportation mode l (just for trucks)
CMLO_ip	Lower bound of dedicated foreign vaccine capacity of p type vaccine on foreign vaccine manufacturing center i
CILO_v	Lower bound of dedicated internal vaccine capacity on internal vaccine manufacturing center v
CPLO_jp	Lower bound of dedicated foreign vaccine capacity of p type vaccine on packaging center j
CDLO_kpm	Lower bound of dedicated foreign vaccine capacity of p type vaccine with inventory tech m on distribution center k
CMMX_ip	Upper bound of dedicated foreign vaccine capacity of p type vaccine on foreign vaccine manufacturing center i
CIMX_v	Upper bound of dedicated internal vaccine capacity on internal vaccine manufacturing center v
CPMX_jp	Upper bound of dedicated foreign vaccine capacity of p type vaccine on packaging center j
CDMX_kpm	Upper bound of dedicated foreign vaccine capacity of p type vaccine with inventory tech m on distribution center k
COV_p	Conversion factor of raw material to foreign vaccine type p
PCK_p	Packaging factor of foreign vaccine type p
PDIS_p	Distributing factor of foreign vaccine type p
COI	Conversion factor of raw material to internal vaccine
DHrt	Vaccine demand at provincial medical center r in period t
CTSlt	Capacity of vehicles of transportation mod l for carrying materials in period t (ton)
IMMAX_ip	Maximum inventory capacity of foreign vaccine type p at foreign vaccine manufacturing center i
IPMAX_jp	Maximum inventory capacity of foreign vaccine type p at packaging center j
IDMAX_kpm	Maximum inventory capacity of foreign vaccine type p at distribution center k with inventory tech m
IVMAX_v	Maximum inventory capacity of internal vaccine type p at internal vaccine manufacturing center v
Environmental parameters
EETVijpl	Total GHG emission of transportation foreign vaccine type p between foreign vaccine manufacturing center i and packaging center j by mod l, km/tone
EETPjkpl	Total GHG emission of transportation foreign vaccine type p between packaging center j and distribution center k by mod l, km/tone
EETHkrpl	Total GHG emission of transportation foreign vaccine type p between distribution center k and provincial medical center r by mod l, km/tone
EETIvrl	Total GHG emission of transportation internal vaccine between internal vaccine manufacturing center v and provincial medical center r by mod l, km/tone
EEPM_ip	Total GHG emission of manufacturing foreign vaccine type p at foreign vaccine manufacturing center i, (unit)
EEPI_v	Total GHG emission of manufacturing internal vaccine at internal vaccine manufacturing center v, (unit)
EEPG_jp	Total GHG emission of packaging foreign vaccine type p at packaging center j, (unit)
EEFP_j	Total GHG emission of building packaging center j,
EEFD_km	Total GHG emission of building distribution center k with inventory tech m,
EEFI_v	Total GHG emission of building internal vaccine manufacturing center v,
EEIM_ip	Total GHG emission of foreign vaccine inventory with type p at foreign vaccine manufacturing center i foreign vaccine manufacturing center i, (unit)
EEIP_jp	Total GHG emission of foreign vaccine inventory with type p at packaging center j, (unit)
EEID_kpm	Total GHG emission of foreign vaccine inventory with type p and inventory tech m at distribution center k, (unit)
EEII_v	Total GHG emission of internal vaccine inventory at internal vaccine manufacturing center v, (unit)
Social parameters
ESTVijpl	Number of local jobs generated per year due to the unit foreign vaccine with type p shipped from foreign vaccine manufacturing center i to packaging center j via transportation mode l
ESTPjkpl	Number of local jobs generated per year due to the unit foreign vaccine package with type p shipped from packaging center j to distribution center k via transportation mode l
ESTHkrpl	Number of local jobs generated per year due to the unit foreign vaccine package with type p shipped from distribution center k to provincial medical center r via transportation mode l
ESTIvrl	Number of local jobs generated per year due to the unit internal vaccine shipped from internal vaccine manufacturing center v to provincial medical center r via transportation mode l
ESMM_ip	Number of local jobs generated per year due to the unit production of foreign vaccine with type p at foreign vaccine manufacturing center i
ESMI_v	Number of local jobs generated per year due to the unit production of internal vaccine at internal vaccine manufacturing center v
ESPG_jp	Number of local jobs generated per year due to the unit packaging of foreign vaccine with type p at packaging center j
ESPD_kpm	Number of local jobs generated per year due to the unit distributing of foreign vaccine with type p at distribution center k
ESDI_v	Number of local jobs generated per year due to the unit distributing of internal vaccine at internal vaccine manufacturing center v
LSTV_l	Average of lost days per year due to the unit local jobs of vaccine transport process by transportation mod l
LSMM	Average of lost days pear year due to the unit local jobs of vaccine production process
LSPG	Average of lost days per year due to the unit local jobs of vaccine packaging process
LSDS	Average of lost days per year due to the unit local jobs of vaccine distributing process
HSMM_ip	Average fraction of potentially hazardous foreign vaccine type p at foreign vaccine manufacturing center i
HSMI_v	Average fraction of potentially hazardous internal vaccine at internal vaccine manufacturing center v
θ_1	Weighting factor of total number of produced job opportunities
θ_2	Weighting factor of total number of lost days caused from work's damages
θ_3	Weighting factor of total number of potentially hazardous vaccine
Integer variables
RTS_l	Number of purchased vehicle mod l for carrying materials (just trucks)
Binary variables
ZP_j	Equal to 1 if packaging center j is opened; 0 otherwise
ZD_km	Equal to 1 if distribution center k with inventory tech m is opened; 0 otherwise
ZI_v	Equal to 1 if internal manufacture center v is opened; 0 otherwise
Continuous decision variables
CAPMipt	The designed capacity of foreign vaccine type p at foreign vaccine manufacturing center i in time period t
CAPPjpt	The designed capacity of foreign vaccine type p at packaging center j in time period t
CAPDkpmt	The designed capacity of foreign vaccine type p at foreign vaccine manufacturing distribution center k with inventory tech m in time period t
CAPIvt	The designed capacity of internal vaccine at internal vaccine manufacturing center v in time period t
PROMipt	Amount foreign vaccine type p produced at foreign vaccine manufacturing center i in period t
PROIvt	Amount of internal vaccine produced at internal vaccine manufacturing center v in period t
PRODkpmt	Amount of foreign vaccine type p pre-distributed at distribution center k with inventory tech m in period t
PGOPjpt	Amount of foreign vaccine type p packaged at packaging center j in period t
SUPNipt	Amount of raw material foreign vaccine type p purchased at foreign vaccine manufacturing center i in period t
SUPIvt	Amount of raw material internal vaccine purchased at internal vaccine manufacturing center v in period t
SVijplt	The amount of foreign vaccine type p transported from foreign vaccine manufacturing center i to packaging center j with transportation mod l at time phase t
SPjkplt	The amount of foreign vaccine type p transported from packaging center j to distribution center k with transportation mod l at time phase t
SHkrplt	The amount of foreign vaccine type p transported from distribution center k to provincial medical center r with transportation mod l at time phase t
SIvrlt	The amount of internal vaccine transported from internal vaccine manufacturing center v to provincial medical center r with transportation mod l at time phase t
INVMipt	The inventory of foreign vaccine type p at foreign vaccine manufacturing center i in period t
INVPjpt	The inventory of foreign vaccine type p at packaging center j in period t
INVDkpmt	The inventory of foreign vaccine type p at distribution center k with inventory tech m in period t
INVIvt	The inventory of internal vaccine at internal vaccine manufacturing center v in period t
CAPEMipt	The capacity expansion of foreign vaccine type p at foreign vaccine manufacturing center i in period t
CAPEPjpt	The capacity expansion of foreign vaccine type p at packaging center j in period t
CAPEDkpmt	The capacity expansion of foreign vaccine type p at distribution center k with inventory tech m in period t
CAPEIvt	The capacity expansion of internal vaccine at internal vaccine manufacturing center v in period t
SQMipt	The shortage of foreign vaccine type p at foreign vaccine manufacturing center i in period t
SQPjpt	The shortage of foreign vaccine type p at packaging center j in period t
SQDkpmt	The shortage of foreign vaccine type p at distribution center k with inventory tech m in period t
SQIvt	The shortage of internal vaccine at internal vaccine manufacturing center v in period t
COSFIX	Total fixed opening costs
COSVAR^t	Total variable opening costs in periods t
COSPRO^t	Total production costs in periods t
COSPKG^t	Total packaging costs in periods t
COSSU^t	Total raw material supplying costs in periods t
COSDIS^t	Total distribution costs in periods t
COSINV^t	Total inventory holding costs in periods t
COSTR^t	Total material transportation costs in periods t
COSSHG^t	Total shortage costs in periods t
COSEXP^t	Total capacity expansion costs in period t
COSINT^t	Total purchasing road equipment costs in period t
CF^t	Total cost in period t
EFTR^t	Total amount of GHG emission impact for materials transportation in period t
EFIN^t	Total amount of GHG emission impact for vaccine inventory in period t
EFPR^t	Total amount of GHG emission impact for vaccine production in period t
EFPG^t	Total amount of GHG emission impact for packaging in period t
EFFR	Total amount of GHG emission impact for opening facilities
NUMW	Total number of local jobs generated
TRSC^t	Total number of local jobs generated in vaccine transportation section in period t
MNSC^t	Total number of local jobs generated in vaccine production section in period t
PGSC^t	Total number of local jobs generated in vaccine packaging section in period t
DSSC^t	Total number of local jobs generated in vaccine distributing section in period t
LOSC	Total average of lost days in total local jobs
HZSC	Total average fraction of potentially hazardous in vaccine production