key: cord-0696349-346aosj8
authors: Yin, Xuecheng; Büyüktahtakın, İ. Esra; Patel, Bhumi P.
title: COVID-19: Data-Driven Optimal Allocation of Ventilator Supply under Uncertainty and Risk
date: 2021-12-01
journal: Eur J Oper Res
DOI: 10.1016/j.ejor.2021.11.052
sha: cc455ddaab773df859f128f6aecb1080a4a86f98
doc_id: 696349
cord_uid: 346aosj8

This study presents a new risk-averse multi-stage stochastic epidemics-ventilator-logistics compartmental model to address the resource allocation challenges of mitigating COVID-19. This epidemiological logistics model involves the uncertainty of untested asymptomatic infections and incorporates short-term human migration. Disease transmission is also forecasted through a new formulation of transmission rates that evolve over space and time with respect to various non-pharmaceutical interventions, such as wearing masks, social distancing, and lockdown. The proposed multi-stage stochastic model overviews different scenarios on the number of asymptomatic individuals while optimizing the distribution of resources, such as ventilators, to minimize the total expected number of newly infected and deceased people. The Conditional Value at Risk (CVaR) is also incorporated into the multi-stage mean-risk model to allow for a trade-off between the weighted expected loss due to the outbreak and the expected risks associated with experiencing disastrous pandemic scenarios. We apply our multi-stage mean-risk epidemics-ventilator-logistics model to the case of controlling COVID-19 in highly-impacted counties of New York and New Jersey. We calibrate, validate, and test our model using actual infection, population, and migration data. We also define a new region-based sub-problem and bounds on the problem and then show their computational benefits in terms of the optimality and relaxation gaps. The computational results indicate that short-term migration influences the transmission of the disease significantly. The optimal number of ventilators allocated to each region depends on various factors, including the number of initial infections, disease transmission rates, initial ICU capacity, the population of a geographical location, and the availability of ventilator supply. Our data-driven modeling framework can be adapted to study the disease transmission dynamics and logistics of other similar epidemics and pandemics.

The world is undergoing a major health crisis, which has now turned into a pandemic. The Coronavirus Disease 2019 , first detected in Wuhan city of China at the end of 2019, has had a devastating effect on human life and economies in all parts of the world. Countries worldwide enforce lockdown and quarantine rules to slow down the spread of the virus. For instance, a 76day strict lockdown was imposed starting on January 23, 2020, public transport was suspended, school re-openings after the winter vacation were delayed, and population movements were severely curtailed in Wuhan. Due to the quick response and successful interventions, China and some other Asia-Pacific countries have managed to control the pandemic rapidly and effectively (Burki, 2020) . New Zealand imposed a five-week nationwide lockdown before returning to close-to-normal. The country closed its borders to foreigners for more than a year, reducing their cases and deaths considerably. On the contrary, a slow reaction delaying interventions, such as mask and social distancing, have aggravated the COVID-19 death toll in some other world countries, such as the United States and the United Kingdom (UK) (Lilleker, 2021; CNN, 2021) . South Korea and Ghana gave consistent messages highlighting the risks of the new pandemic and how it could be mitigated. However, nations such as Brazil, the UK, and India gave out inconsistent messages about the threat of COVID-19 and encouraged complacency less well, worsening the pandemic's impact (CNN, 2021) .

The lockdown, imposing travel restrictions, and social distancing have also severely affected the economy, from small-scale industries to stock prices and international trading. The virus has a high transmission rate, causing more than 104.7 million cases globally, of which 2.3 million people have died by mid-February 2021 (JHU, 2020) . The continuous increase seen in coronavirus cases has made a worldwide scarcity of essential resources, such as ventilators, Intensive Care Unit (ICU) beds, Personal Protective Equipment (PPE), and masks. Effective, sufficient, and timely delivery of those critical resources to serve COVID-19 patients has been a major challenge faced by nations worldwide during the pandemic.

COVID-19 is primarily an acute respiratory disease. Ventilator incubation delivers high oxygen concentrations while removing carbon dioxide and reduces the risk of hypoxia for COVID-19 patients . The standard Acute Respiratory Distress Syndrome protocol mandates that the most severe COVID-19 patients, who constitute 5% of all COVID-19 patients, should receive ventilator support (Bein et al., 2016) . As a result, the life of many COVID-19 patients depends on the use of ventilators. The shortage of supplies and uncertainty in disease transmission has affected the proper allocation of ventilators, causing immense distress to the healthcare system. Due to ventilator shortages worldwide during the pandemic's peak times, hospital officials have had to make life-altering resource allocation decisions and prioritized the care of COVID-19 patients (Ranney et al., 2020) . To tackle ventilator shortages and reduce the number of COVID-19-related deaths, studies have come up with new approaches for ventilator distribution. For example, Ranney et al. (2020) suggest that the demand for ventilators can be fulfilled by the government by allowing other industries to come together and help medical industries to cater to the needs of the ventilators. Another study by Castro et al. (2020) suggests that the government in Brazil should start thinking about expanding the resource capacity rather than only focusing on the allocation of the available resources for controlling COVID-19. White and Lo (2020) develop a framework for the distribution of ICU beds and ventilators depending on the priority scores using a scale of 1 to 8 based on patients' likelihood of survival and ethical considerations.

Operations Research (OR) methods have been widely used to determine optimal resource allocation strategies to control an epidemic or pandemic. Several studies have used multi-period OR models to optimize the allocation and redistribution of ventilators (see, e.g., Mehrotra et al. (2020) , Bertsimas et al. (2020) , and Blanco et al. (2020) ). Other OR research models that study the epidemic diseases and resource allocation mainly focus on the logistics and operation management to control the disease in optimal ways Zaric and Brandeau, 2001; Yin and Büyüktahtakın, 2021a; Kaplan et al., 2003; Tanner et al., 2008; Coşgun and Büyüktahtakın, 2018) . We refer the reader to excellent reviews of Dasaklis et al. (2012) and Queiroz et al. (2020) for a discussion of OR models for epidemic resource allocation.

While OR has been an extremely useful tool for effective and timely allocation of resources as a response to epidemics, none of the former work has considered the ventilator allocation problem using a risk-averse spatio-temporal stochastic programming model under uncertainty of asymptomatic infections. People move between regions, states, and countries, which aggravates the disease transmission to the other areas. Evaluating undetected or asymptomatic individuals is critical for determining disease dynamics because asymptomatic individuals move around and unknowingly infect other individuals (McCrimmon, 2021) . Thus, the short-term migration of people is a critical factor that needs to be considered to forecast the transmission of COVID-19 realistically. However, the short-term migration rate is hard to predict and is affected by interventions and human behaviors. Furthermore, disease transmission rates are not constant and rather evolve over time with government interventions, such as the lockdown or social distancing measures. This change in the transmission rates also should be considered in a realistic model. To our knowledge, none of the former OR ventilator allocation models have integrated the epidemiological aspects of the disease and resource allocation challenges in one optimization model.

In this paper, we address the limitation of realistically forecasting the transmission of COVID-19 and build a risk-averse multi-stage stochastic epidemics-ventilator-logistics programming model to study the ventilator allocation for the treatment of severe COVID-19 patients. Our model considers the uncertainty of untested asymptomatic individuals during the transmission of COVID-19 and involves various pandemic scenarios during each time stage of the planning horizon. Our model also incorporates the short-term migration between the highly-impacted regions while using changing transmission rates under various non-pharmaceutical intervention measures. The model optimizes the distribution of ventilators while minimizing the total expected number of infected and deceased people. We calibrate, validate, and test our epidemiological ventilator allocation model using COVID-19 data collected during the pandemic's early stages.

Former research presents bounds and approximation methods for risk-averse multi-stage stochastic mixed-integer programming problems, see, e.g., Sandıkçı andÖzaltın (2017), Mahmutogulları et al. (2018) , and Bakir et al. (2020) . Büyüktahtakın (2021) presents a new scenario sub-problem and defines stage-t scenario dominance concept for the first time. The author also defines lower and upper bounds on the problem based on scenario sub-problems and presents scenario dominance decomposition cuts based on a partial ordering of scenarios. Bushaj et al. (2021a) later extend the scenario dominance definition and decomposition cuts of Büyüktahtakın (2021) under the case of endogenous uncertainty and show their computational benefit for solving a forest insect infestation surveillance and control problem. Inspired by the approach of Büyüktahtakın (2021), we define region-based and risk-based sub-problems and new lower and upper bounds using those sub-problems. Computational results also show that the proposed bounds significantly reduce the optimality and relaxation gaps.

This paper is organized into six sections. Section 2 presents a review of the articles that estimate the transmission of COVID-19, evaluate various interventions used to control the disease, and optimize resource allocation during an epidemic or pandemic, and summarizes our paper's contributions. Section 3 presents the model notations, the compartmental sub-model, time and space-varying transmission sub-model, the description of uncertainty and scenario tree, the risk measure used in the model, and the multi-stage stochastic programming model formulation that integrates all elements described in this section. Section 4 presents a new region sub-problem and lower and upper bounds to tackle the problem difficulty. Section 5 provides the data used in the model, including population and migration data, epidemiological data, initial infection data, and capacity and cost data. Section 6 presents the results of the model's application to the case of COVID-19 control in the counties of New York and New Jersey states. Section 7 discusses key insights from the results and provides future research directions.

This section presents a review of the articles that study mathematical models to estimate the transmission rate of COVID-19, evaluate various interventions used to control the disease, and optimize resource allocation during the pandemic. The literature covers various methods, including compartmental models, such as the Susceptible (S)-Exposed (E)-Infectious (I)-Recovered (R) [SEIR] formulations, simulations, optimization models, stochastic, statistical, or probabilistic approaches, and network or graph-theory models, and is summarized in Table 1 .

Compartmental, Simulation and Network Models. Several studies on COVID-19 modeling use SEIR-type models with simulations, such as the Monte Carlo Simulation, to predict the disease transmission rate and analyze the effectiveness of proposed non-pharmaceutical interventions and resource allocation strategies (see, e.g., Chatterjee et al. (2020) , Weissman et al. (2020) , Wang et al. (2020) , and Li et al. (2020) ). For example, Li et al. (2020) develop a novel epidemiological model called DELPHI based on an SEIR model to estimate the effectiveness of government interventions and the effects of under-detection of confirmed cases. They find that the world-widely implemented travel restriction policies, social distancing, and mass gathering restrictions play a crucial role in reducing the infection rate. Ku et al. (2020) develop an empirical approach based on a Bass Susceptible-Infected-Recovered model to predict specific transmission parameters, exogenous forces of infection, and effective population sizes to determine the reproductive number (R 0 ) for COVID-19 transmission in Chinese provinces.

Stochastic compartmental models have also been used to analyze the uncertain transmission dynamics of COVID-19. For example, Kucharski et al. (2020) develop the stochastic dynamic Susceptible, Isolated, Infected, Recovered, and Exposed model to derive the transmission rate using COVID-19 cases in Wuhan as well as the associated international infections. The initial reproductive number of 2.35 was reduced to 1.05 after imposing travel restrictions in Wuhan. Kretzschmar et al. (2020) develop a stochastic mathematical model to understand the impact of time delays in testing and isolation on the reproductive number of COVID-19 transmission. Differential equations have been mainly used to solve compartmental disease models (see, e.g., Ambikapathy and Krishnamurthy (2020) , Wei et al. (2020) , Zeb et al. (2020) , Tuan et al. (2020), and Wang et al. (2020) ). For example, Roberts et al. (2020) come up with a series of differential equations that are designed to compare the two scenarios based on delaying of ICU bed shortageeffects of hospitalizing fewer COVID-19 patients versus increasing the ICU bed capacity.

The literature also includes many studies that present statistical approaches, such as logistic regression, generalized additive, and time-series models to predict the growth of the disease over time and implement resource allocation strategies accordingly during the pandemic (see, e.g., Murray (2020), Manca et al. (2020) , and Katris (2020)). Graph-theoretical methods have also been quite useful in studying the COVID-19 transmission patterns and finding out the practices to reduce the transmission rate (see, e.g., Hu et al. (2020 ), Loeffler-Wirth et al. (2020 , and Badr et al. (2020) ). For example, Badr et al. (2020) model the transmission rate of COVID-19 based on mobile phone mobility and support the role of social distancing in combating COVID-19.

Optimization Models. Optimization models have also been widely studied for resource allocation in the fight against COVID-19. Queiroz et al. (2020) provide a systematic review of various supply chain and logistics approaches for optimizing the distribution of critical resources amid COVID-19. To tackle the shortage of ventilators, Mehrotra et al. (2020) develop a two-stage stochastic programming model, optimizing ventilator allocation during the pandemic under various demand scenarios. The authors find that when 60% of the ventilator inventory is allocated to non-COVID-19 patients, there is no shortfall. In comparison, when 75% of the stock is allocated to the non-COVID-19 patients, a shortfall in the supply of the ventilators to COVID-19 patients occurs. Also, they find that it is essential to ramp up the production of the ventilators to meet the additional requirements of the ventilators that might come up during the peak times of the pandemic. Lacasa et al. (2020) come up with an algorithm for optimizing the allocation of the ventilators and ICU beds and validate their algorithm during the peak and declining times of the pandemic based on the data from the United Kingdom and Spain. Bertsimas et al. (2020) develop a four-step approach, combining descriptive, predictive, and prescriptive analytics and propose an optimization model for the re-allocation of the ventilators throughout the U.S. during the COVID-19 pandemic. Blanco et al. (2020) present a two-stage stochastic mixed-integer programming model, which minimizes the expected non-covered demand, using robust objective functions of type min-max regret. Billingham et al. (2020) present a network optimization model to tackle the problem of scarce ventilator distribution. Parker et al. (2020) develop mixed-integer programming and robust optimization models to redistribute patients instead of resources, such as ventilators among different hospitals under demand uncertainty. Govindan et al. (2020) develop a practical decision support system hinging on the knowledge of the physicians and the fuzzy interference system (FIS) to help manage the demands of essential hospital services in a healthcare supply chain, to break down the pandemic propagation chain, and reduce stress among the health care workers. Huang et al. (2017) introduce a two-stage forecasting and optimization framework for optimizing stockpiles of mechanical ventilators, which are critical for treating hospitalized influenza patients in respiratory failure under a pandemic situation. They also incorporate their model into a web-based decision-support tool for pandemic preparedness and response. Yin and Büyüktahtakın (2021b) develop a multi-stage stochastic programming vaccine allocation model and apply it to the case of the Ebola virus disease in the Democratic Republic of the Congo. The authors also define two new equity metrics for fair resource allocation. Bushaj et al. (2021b) present an integrated simulation-deep reinforcement learning (SiRL) framework to determine best intervention strategies over time for epidemic decision-making, while Yin et al. (2021) present an agent-based simulation-optimization approach to solve the vaccine center location vaccine allocation problem with an application to controlling COVID-19.

In summary, former stochastic programming approaches on ventilator allocation in a pandemic situation have involved a time domain of only two stages and have not integrated an epidemic model within the stochastic program. Furthermore, the mathematical models on the forecast of COVID-19 do not include the uncertainty of untested asymptomatic infections. They do not incorporate the impact of short-term migrations on COVID-19 transmission in an epidemiological model. Also, former studies on COVID-19 modeling and logistics have omitted the time consistency of the risk for making decisions over multiple stages of a stochastic program under extreme pandemic scenarios. Our methodological, modeling, and applied contributions to the epidemiology and OR literature are summarized below.

Methodological contributions: We define region-and risk-based sub-problems. Then using those sub-problem solutions, we generate new lower and upper bounds on the original problem optimal objective function value. We also present proofs for those results and bounds. Our computations demonstrate that the proposed region-based upper and lower bounds significantly reduce the optimality and relaxation gaps, respectively. We also compare and discuss the efficiency of region sub-problem bounds with respect to that of the scenario sub-problem bounds of Büyüktahtakın (2021) .

Modeling contributions: First, to our knowledge, this is the first study that addresses the optimal distribution of ventilators to control a pandemic in a multi-stage stochastic mean-Conditional Value at Risk (CVaR) model. This model includes many realistic effects critical in the COVID-19 pandemic, including untested asymptomatic infections, human movement among multiple regions, and evolving transmission rates under non-pharmaceutical intervention measures. Considering multiple stages in a multi-stage SP is essential to capture uncertain disease dynamics over multiple time periods since disease transmission dynamically changes over time. Second, we consider the uncertainty of the proportion of untested asymptomatic infections at each stage and integrate this unknown dimension of the pandemic by generating a multi-stage scenario tree. As opposed to robust optimization, which focuses on the worst-case outcomes assuming a hostile environment, in stochastic programming, decision makers consider all possible outcomes (e.g., scenarios) by assigning weights to outcomes and optimizing some aggregated measure of performance over all likely outcomes (Defourny et al., 2012) . Thus, our multi-stage stochastic programming approach to modeling the uncertainty in the proportion of asymptomatic infections is more advantageous and provides less conservative results compared to robust optimization (Yin and Büyüktahtakın, 2021a) .

Third, we present a new susceptible (S)-tested infected (I)-untested asymptomatic (X)-hospitalized (H)-ICU (C)-recovered (R)-death (D) compartmental disease model specialized for COVID-19, and also incorporate short-term human migration among multiple regions into this epidemiological model. Fourth, we derive a new time-and space-varying disease transmission formulation, which takes into account the impact of government interventions on transmission rates. Fifth, we formulate a budget-constrained ventilator allocation logistics model. Sixth, we incorporate a time-consistent CVaR risk-measure and the expectation criterion in the objective function to alleviate the impacts of extreme pandemic scenarios. Lastly, we integrate all those elements into one epidemics-ventilator-logistics mathematical formulation, which minimizes the number of infections and deceased individuals under different intervention strategies while determining the optimal timing and location of resources (ventilators) allocated. Our model combines the forecast of the transmission of COVID-19 and the determination of optimal ventilator allocation strategies in one formulation. Accordingly, the decision-maker can evaluate possible outcomes of wait-and-see decisions while foreseeing how the disease could progress in each time period.

Applied Contributions: We apply our general multi-stage mean-risk epidemics-ventilator-logistics model to the case of controlling COVID-19 in highly-impacted counties of New York and New Jersey. We collect real data from various resources and provide researchers with compact epidemiological, population, and logistics-capacity data for COVID-19. Using this data, we calibrate, validate, and test our model, which could be used as a decision support tool for fighting against COVID-19. Our model can also be adapted to study other similar diseases' transmission dynamics and logistics.

This study provides optimal risk-averse ventilator allocation policies under different risk levels that the decision-maker can apply to control COVID-19. Based on our results, we offer the following recommendations to inform resource allocation policies under a pandemic:

The short-term movement of people influences the number of new infections even if the disease transmission rate stays the same. The number of treated people in the ICU may stay at the capacity limit under different intervention strategies because this value depends on the minimum number of patients who require a ventilator for treatment and the scarce ventilator supply. There is also a lag time to observe the impacts of government non-pharmaceutical interventions on the number of hospitalized and deceased individuals, as well as ICU patients. "Lockdown" is the best strategy to control COVID-19. However, the "Mask and Social Distance" intervention following a certain period of "Lockdown" is the second-best choice, considering the need for opening facilities and businesses. The ventilator allocation under different budget levels and scenarios indicates that the number of ventilators allocated to each region depends on various factors, such as the number of initial infections, initial disease transmission rates, initial ICU capacity, and the population of a geographical location. Overall, the region with a high initial transmission rate and low initial ICU capacity has more ventilators allocated. Under a limited budget level, the region with low initial infections and low initial ICU capacity gets more ventilators allocated when the transmission scenario changes from "All Low" to "All Medium" and "All High." Independent from the transmission scenarios, the region with a high initial transmission rate has more ventilators allocated when the budget level increases. There are also increasing trends in the number of ventilators allocated when either the budget level increases or the cost for each ventilator reduces. A large-enough budget also provides some flexibility in delaying ventilator allocation to some regions. In contrast, all of the ventilators are allocated at the first two stages under a limited budget. Considering risk in decision-making improves the confidence level for reducing the loss of lives under risky pandemic scenarios. However, a risk-averse decision-maker should also expect a possible increase in the number of infections and deaths while mitigating disastrous scenarios.

This section discusses important features and submodels used in the risk-averse multi-stage stochastic epidemics-ventilator-logistics compartmental model and then presents its formulation (5a)-(5w). Specifically, subsection 3.1 presents the description of the notations used in the model, subsection 3.2 describes compartment flow of population dynamics in the mathematical model, subsection 3.3 introduces the formulation for time and space-varying transmission rates in the model, subsection 3.4 explains the uncertainty and scenario tree generation scheme, subsection 3.5 states the specific features and assumptions made in the mathematical model, and subsection 3.6 provides a brief description of the measure of multi-stage risk used in the model, CVaR. Then subsection 3.7 presents the mathematical formulation for our epidemics-ventilator-logistics model that integrates all elements described in this section.

Below we provide notations used for the rest of the paper. r: Index for region where r ∈ R.

ω: Index for scenario, where ω ∈ Ω.

State Variables:

S ω j,r : Susceptible individuals in region r at stage j under scenario ω.

I ω j,r : Tested symptomatic infected individuals in region r at stage j under scenario ω.

X ω j,r : Untested asymptomatic infected individuals in region r at stage j under scenario ω.

H ω j,r : Hospitalized individuals in region r at stage j under scenario ω.

C ω j,r : Individuals treated in the intensive care unit (ICU) in region r at stage j under scenario ω.

R ω j,r : Recovered individuals in region r at stage j under scenario ω.

F ω j,r : Deceased individuals in region r at stage j under scenario ω.

O ω j,r : Number of tested symptomatic infected individuals admitted to the hospital in region r at stage j under scenario ω.

I ω j,r : Number of tested symptomatic infected individuals who cannot be admitted to the hospital due to limited capacity in region r at stage j under scenario ω.

C ω j,r : Number of individuals admitted to ICU in region r at stage j under scenario ω.

K ω j,r : Number of hospitalized individuals not admitted to the ICU due to the limited availability of ventilators in region r at stage j under scenario ω.

U ω j,r : Number of cumulative ventilators (ICU capacity) in region r at stage j under scenario ω.

I ω j,r : Number of infections caused by short-term migration in region r at stage j under scenario ω.

λ 1,r : Recovery rate of tested symptomatic infected individuals in region r.

λ 2,r : The death rate of tested symptomatic infected individuals in region r. λ 3,r : Hospitalization requirement rate of tested symptomatic infected individuals in region r.

λ 4,r : Recovery rate of the hospitalized individuals in region r.

λ 5,r : Death rate of hospitalized individuals in region r. λ 6,r : Ventilator requirement rate of hospitalized individuals in region r.

λ 7,r : Recovery rate of ICU patients in region r.

λ 8,r : Death rate of ICU patients in region r.

λ 9,r : Recovery rate of untested asymptomatic individuals in region r.

σ 1,j,r : Transmission rate of tested symptomatic infected individuals in region r at stage j.

σ ω 2,j,r : Proportion of untested asymptomatic infections in region r at stage j under scenario ω.

T ω j,r : Hospital capacity in region r at stage j under scenario ω. π r : Initial number of susceptible individuals in region r. ϖ r : Initial number of tested infected individuals in region r. φ r : Initial number of untested asymptomatic infected individuals in region r. κ r : Initial number of hospitalized individuals in region r. υ r : Initial number of ICU patients in region r. ϑ r : Initial number of recovered individuals in region r. τ r : Initial number of deceased individuals in region r.

α: Confidence level of value-at-risk, where α ∈ [0, 1). ϕ: Non-negative risk preference parameter or mean-risk trade-off coefficient.

η ω j : Value at risk for each stage j under scenario ω.

z ω j : Value exceeding the value-at-risk at the confidence level α at stage j under scenario ω.

n: The serial number of nodes in the scenario tree, where n ∈ N . t(n): The corresponding stage that node n marked in the scenario tree.

The set of scenarios that pass through node n.

y ω j,r : Number of ventilators allocated to region r at the end of stage j under scenario ω.

Q ω j,r : Auxiliary variable to be substituted with λ 3,r I ω j,r .

W ω j,r : Auxiliary variable to be substituted with T ω j,r − H ω j,r .

V ω j,r : Auxiliary variable to be substituted with λ 6,r H ω j,r .

G ω j,r : Auxiliary variable to be substituted with U ω j,r − C ω j,r .

m ω j,r : Binary variable, which takes the value 0 if the number of infected individuals who need to be hospitalized is restricted by the number of available beds in hospitals, and the value 1 if the number of beds in hospitals is sufficiently large to hospitalize all infected individuals who need to be hospitalized.

d ω j,r : Binary variable, which takes the value 0 if the number of hospitalized individuals need to be treated in ICUs is restricted by the number of available ventilators, and the value 1 if the number of ventilators is sufficiently large to treat all hospitalized individuals who need to be treated in ICUs.

A LB : Lower bound of λ 6,r H ω j,r .

B U B : Upper bound of λ 3,r I ω j,r .

B LB : Lower bound of λ 3,r I ω j,r .

This subsection describes the compartment flow of population dynamics at each stage in each region. Figure 1 shows the transmission dynamics of COVID-19 in each region r at each time period j for a particular scenario ω. In this figure, susceptible individuals (S) can be infected and become infected (either symptomatic or asymptomatic). Asymptomatic infections (X) may have slight or no symptoms throughout the infection period and will recover with a rate of λ 9,r . Tested symptomatic infections (I) may recover or die with rates of λ 1,r and λ 2,r , respectively, if they are not treated in the hospital (H). Note thatÎ ω j,r with an incoming dashed arc to the I compartment represents the number of infected people coming into the region r at stage j under scenario ω from neighboring regions. Tested infected individuals (I) move to the hospital (H) compartment, depending on the number of tested infections (I) and available hospital capacity (we use λ 3,r to represent the associated non-constant transmission rate). Some of the treated infected people in the hospital (H) will recover with a rate of λ 4,r , where R represents recovered individuals. The Figure 1 : One-step COVID-19 compartmental model. Note that λ 3,r and λ 6,r are not constant transition rates; and they depend on the available hospital and ICU capacity, respectively.

situations of some patients in the hospital (H) may worsen, and thus they may be transferred to the intensive care unit (we use C to represent ICU), and those individuals need ventilators for their treatment (similar to λ 3,r , λ 6,r represents the non-constant rate from (H) to (C) since this number depends on the minimum number of patients who need to be treated in ICU and the available ventilator supply in the ICU).

Similar to the case of admittance into the hospital, the number of hospitalized patients transferred to an ICU at each time period is equal to the minimum of the number of patients who need to be transferred to an ICU and the number of available ventilators. The patients who are not able to receive the treatment in the ICU due to the limitation on the number of available ventilators may die at a rate of λ 5,r , where F represents deceased individuals. After being treated in the ICU, some of the patients may recover with a rate of λ 7,r , while others may die with a rate of λ 8,r . Different from a typical compartmental model, the transfer rate from I to H and H to C is not a constant, and it depends on the available capacity in the H and C compartments, respectively, as discussed above.

This subsection introduces how we model COVID-19 transmission under government interventions in the mathematical model. We formulate the transmission rate σ 1,j,r as a time-and spacevarying parameter, which depends on the government interventions taken at time j and region r. Since the onset of COVID-19, many governments have imposed different intervention strategies to reduce the transmission rate. At a certain stage, each intervention has a different impact on the transmission rate for the next stage. In this paper, we incorporate three main non-pharmaceutical interventions to formulate the time-varying transmission rate -none, mask and social distancing, and lockdown. Under the mask and social distancing, people are allowed to go to public places, but they have to wear masks and keep a 6-feet distance from each other. The lockdown represents that all the public areas are closed, and people must stay where they are. During the lockdown period, people may go out to buy necessities, but they have to wear masks and keep a 6-feet distance from each other.

Let x 1 j,r , x 2 j,r , and x 3 j,r be binary decision variables that correspond to none (i = 1), mask and social distancing (i = 2), and lockdown (i = 3) interventions, respectively, taken at stage j and region r. If x i j,r takes a value of 1, then intervention i is employed; otherwise, it is not employed at stage j and region r. The transmission rate at stage j + 1 in region r is a function of the transmission rate and specific intervention employed at stage j in the same region, as given in the below equations:

where m i j,t represents the percent change in the transmission rate with respect to the binary decision variable x i j,r for intervention i = 1, 2, 3 taken at stage j in region r. Equation (1) shows that the transmission rate at stage j + 1 is a function of the transmission rate at stage j and the intervention strategy i taken at stage j. Equation (2) indicates that only one intervention measure can be taken at each stage j. Equation (3) describes the binary nature of intervention decisions.

The transmission rate in our model is not equal to the basic reproduction number, R 0 . It shows how many new tested infections will be caused by the symptomatic and asymptomatic infections from the previous stage. Since the number of new asymptomatic infections is uncertain, the number of new infections (both symptomatic and asymptomatic) changes under different scenarios even if the transmission rates at each stage j stay the same.

There is a delay in the impact of the government's interventions on the number of infections and the reaction to the test results is also slow. Therefore, we calculate the transmission rate for the first two stages directly using the real data from JHU (2020), independent of the intervention type. Based on the first two-stage transmission rates, we calculate the transmission rates from stages three to five using the formulation (1)-(3) for each intervention strategy. Also, the values of m i j,t are trained using the real data obtained from JHU (2020). As an example, the initial transmission rates for the first two stages in New York and New Jersey and the impacts of government intervention strategies m i j,t are shown in Table 7 under Section 5.2.

In this subsection, we explain how we address the uncertainty of disease transmission in the mathematical model. Data regarding undetected or untested asymptomatic cases is lacking and uncertain. Therefore, we model the uncertainty regarding the proportion of untested asymptomatic infections (σ ω 2,j,r ) by generating a set of scenarios ω ∈ Ω, each representing a specific realization of the uncertain proportion of untested asymptomatic individuals over multiple time periods. Our scenario generation approach is similar to Alonso-Ayuso et al. (2018)'s method developed to model the demand uncertainty in forestry management. Each scenario has a probability of p ω and ω∈Ω p ω = 1. Since data is not available to describe the probability distribution of the uncertain variable (σ ω 2,j,r ), we assume that the uncertain parameter follows a normal distribution. The lower and upper bounds for the proportion of asymptomatic infections are obtained from the study of Meller (2020) . The lower bound value for the random variable is considered as the value of 0.001-quantile, and the upper bound is considered as the value of 0.999-quantile of the normal distribution.

As an example, Figure 2 shows a particular scenario tree for the proportion of untested asymptomatic infections (σ ω 2,j,r ) for a two-stage problem. We consider three realizations at each node of the scenario tree by dividing its normal distribution into three discrete parts [low (L), medium (M), high (H)]. The low and high realizations have a probability of 0.3, and the medium realization has a probability of 0.4. Each path from the root node to the leaf node of the scenario tree represents a scenario ω. The probability of a scenario ω, p ω , is calculated as the multiplication of probabilities on the path for scenario ω. For two stages, 9 (3 2 ) scenarios will be generated in this instance. The non-anticipativity constraints indicate that two scenarios are inseparable at stage j if they share the same scenario path up to that stage. This means that the corresponding decision made at this stage for those two scenarios should also be the same.

The value of the proportion of asymptomatic infections has a mean µ j r and standard deviation σ j r at stage j. We use Q h to represent the value of h-quantile in the normal distribution. For each node n in the scenario tree, the mean value of the low realization is the value of 0.15-quantile (E(µ n r,low |Q 0.001 ≤ µ n r,low ≤ Q 0.30 ) = Q 0.15 ), the mean value of medium realization is the value of 0.50-quantile (E(µ n r,medium |Q 0.30 ≤ µ n r,medium ≤ Q 0.70 ) = Q 0.50 ), and the mean value of high realization is the value of 0.85-quantile (E(µ n r,high |Q 0.70 ≤ µ n r,high ≤ Q 0.999 ) = Q 0.85 ). For node 0 in our example, the proportion of untested asymptomatic infections at stage j = 0 has µ 0 r = 0.26 and σ 0 r = 0.05. The low, medium, and high realizations at node 0 in stage j = 0 and nodes 1 and 3 in stage j = 1 are given in Table 2 below. According to the distributions presented in Table 2 , the proportion of untested asymptomatic infections in stage 1 is realized as 0.21 (Low) at node 1, 0.26 (Medium) at node 2, and 0.31 (High) at node 3.

Decisions at stage j ∈ J = {0, ...,J} are made based on the available information up to stage j, which is also known as the non-anticipativity requirement. In our multi-stage stochastic program, we have two types of decision variables, namely here-and-now and wait-and-see variables. Hereand-now variables refer to ventilator allocation decision variables (y ω j,r ) that must be decided at the beginning of each time stage j prior to the realization of the uncertain parameter (σ ω 2,j,r ) at that stage. Wait-and-see variables, on the other hand, refer to State Variables presented under Section 2.1 (e.g., S, I, and X) that are to be determined at the end of each time stage j after ventilator allocation decisions at that stage are made and the outcome of the uncertain parameters are known. For example, at the beginning of stage j we take a ventilator allocation action before the uncertain parameter σ ω 2,j,r is known at stage j. Once the uncertain parameter value is realized at stage j, we update the state variables at stage j + 1 based on the optimal treatment actions and the ventilator allocation decisions made up to stage j, and then take another ventilator allocation action before the uncertainty in stage j + 1 is realized. This is followed up by updating the value of state variables based on the optimal treatment actions and the ventilator allocation decisions made up to that point in the subsequent stages.

The decision process in a multi-stage stochastic problem is equivalent to sequential decision making over multiple periods where here-and-now and wait-and-see decisions follow each other, as illustrated in an example in Figure 3 . In this example figure, we represent here-and-now and wait-and-see decisions to be taken at each stage j = {1, ...,J − 1} as y ω j,r and I ω j,r , respectively, and the uncertain parameter in stage j and in region r is represented by ξ j,r . In each time stage j, the realization of ξ j,r is known only after the decisions at stage j, y ω j,r , have already been made. Once ξ j,r is realized, we make the state decisions represented by I ω j+1,r at stage j + 1. Figure 2 and the associated node of the uncertain parameter realization.

Low (realized node) Medium (realized node) High (realized node)

Node 0 Distribution 0.21 (node 1) 0.26 (node 2) 0.31 (node 3)

Node 1 Distribution 0.17 (node 4) 0.21 (node 5) 0.25 (node 6)

Node 3 Distribution 0.12 (node 10) 0.31 (node 11) 0.50 (node 12)

Since the transmission of COVID-19 is affected by many factors, data to calibrate some of the model parameters, such as the impact of human mobility, is either lacking or inaccurate. Therefore, we incorporate some important features and make some assumptions in the model formulation.

Important features. First, we consider the impact of different intervention strategies on the disease transmission rate and adjust the short-term migration population depending on the intervention strategy. For instance, under the lockdown strategy, we assume that the short-term migration among each county is zero. Under mask and social distancing strategies, the short-term migration population among each county is reduced to 60% of the original value, as estimated from the study of Lee et al. (2020) . Second, we incorporate the cost for purchasing ventilators to provide a capacity limitation on the total number of ventilators that could be allocated for treating COVID-19 patients. Since there are significant fluctuations in the ventilator prices (Glass, 2020), we consider the minimum purchase price for each ventilator acquired. Third, we train the real data to determine the impact rate of each intervention strategy on the disease transmission rate. The trained value of the impact of interventions can only be used in the regions considered in our case study since all the selected counties in New York and New Jersey are geographically close to each other, and thus interventions have similar social effects. However, for example, the impact of intervention strategies in a rural area may be different from those taken in a city. To estimate the transmission of COVID-19 in other U.S. regions, the model should be re-trained using the associated data.

Assumptions. First, the proportion of untested asymptomatic people is an uncertain parameter in our model. Studies regarding the asymptomatic infections specific to New York City region are limited. Since studied counties in New Jersey and New York are geographically and socially close to each other, the proportions of untested asymptomatic infections at each stage j under scenario ω are set to be the same for each region r. We make this assumption based on the characteristics of these regions. Many people who live in these counties work in New York City, and the short-term migration among each region is quite frequent. Thus, we use the same lower and upper bound of untested asymptomatic infections for these regions. Although the asymptomatic proportion is the same under each scenario, the number of untested asymptomatic infections is different in each region because of the difference in the transmission rate and the initial number of infections. Second, the model considers allocating newly purchased ventilators for the treatment of COVID-19 patients instead of re-allocating existing ventilators from other counties or states since the demand for ventilators during the disease's peak periods is high for all the counties and states, and there is a lead time for transfer of the ventilators between the states that are far from each other. Here, we also assume a central decision maker is entitled to allocate a given supply of ventilators to multiple regions. Third, the infected individuals who cannot be treated in the hospital (both severe and less severe) due to the limited capacity have the same death rate as people who cannot be treated in the ICU due to the limited ventilator supply. This is because the health condition of both groups may worsen without professional treatment. Fourth, we assume that "I" compartment represents tested infected individuals, including all symptomatic and part of the asymptomatic infections. Thus, all the symptomatic infections and part of the asymptomatic infections are tested. The uncertain parameter is the proportion of the untested asymptomatic infections. It is natural for humans to protect themselves as they learn more about an epidemic due to their survival instincts. Thus, fifth, we assume that people react to the pandemic by anticipating the government's interventions and may start social distancing and quarantining days or weeks before an intervention is imposed, as also reported in multiple studies, see, e.g., Zhang et al. (2020) and Fischer et al. (2020) . Therefore, the transmission rate with either doing nothing or mask and social distancing shows a decreasing trend in later stages of the pandemic due to physical distancing among people. Sixth, we use each county's ICU capacity from JHU (2020) as the initial ventilator availability. We assume that non-COVID-19 patients use 60% of this capacity (Mehrotra et al., 2020) . Thus, only 40% of the initial ICUs are available for treating COVID-19 patients. Lastly, the incubation period of COVID-19 ranges from 2-14 days (CDC, 2021). We assume that the "infected" compartment includes the incubation period. Thus, we use a two-week period for each time stage to make sure people who are infected but still in the incubation period at the beginning of each stage are counted as infected at the same stage.

The α-quantile of the cumulative distribution of a random variable z, inf η {η ∈ R : F z (η) ≥ α}, is defined as the value-at-risk (VaR) at the confidence level α ∈ {0, 1} and denoted by VaR α (z). The conditional expected value that exceeds the VaR at the confidence level α is called conditional value-at-risk (CVaR), defined as CVaR α (z) = E(z | z ≥ VaR α (z)). For a minimization problem, VaR α is the α-quantile of the distribution of the cost, and it provides an upper bound on the cost that is exceeded only with a small probability of 1 − α. CVaR α measures an expectation of the cost that is more than VaR α , and can be calculated as an optimization problem as follows (Rockafellar and Uryasev, 2002) :

where (a) + := max(a, 0) for any a ∈ R. We formulate our model as a mean-risk minimization problem:

where E(f (x, ω)) is the expected cost function over the scenarios ω ∈ Ω, CVaR α represents the conditional value-at-risk at α, and ϕ ∈ [0, 1] is a non-negative weighted risk coefficient, and it can be adjusted for a trade-off between optimizing an expectation value and the level of risk taken.

Time Consistency. Time consistency is considered as a critical issue when modeling a risk-averse multi-stage stochastic program. Time consistency means that if you solve a multi-stage stochastic programming model today, you should get the same solution if you resolve the problem tomorrow, given the information that is observed and decided today. We consider a nested risk measure, expected conditional value-at-risk (E-CVaR), as defined in Homem-de Mello and Pagnoncelli (2016) since it is shown to satisfy the time consistency of multi-stage stochastic programs. The E-CVaR can be linearized and formulated as a linear stochastic programming model. In the following section, we will utilize the E-CVaR as a risk measure in our formulation.

The mathematical formulation for our risk-averse multi-stage stochastic epidemics-ventilatorlogistics model (5a)-(5w) is given below.

Epidemics-Ventilator-Logistics Model Formulation:

C ω 0,r = υ r , R ω 0,r = ϑ r , F ω 0,r = τ r , ∀r ∈ R, ∀ω ∈ Ω,

S ω (j+1),r = S ω j,r − σ 1,r (I ω j,r + X ω j,r ) − σ 1,r (I ω j,r + X ω j,r )

I ω (j+1),r = I ω j,r +Î ω 1,j,r + σ 1,r (I ω j,r + X ω j,r ) − λ 1,r I ω j,r − λ 2,r I

R ω (j+1),r = R ω j,r + λ 1,r I ω j,r + λ 9,r X ω j,r + λ 4,r H ω j,r + λ 7,r C ω j,r j ∈ J \ {J}, r ∈ R, ∀ω ∈ Ω,

F ω (j+1),r = F ω j,r + λ 2,r I ω j,r + λ 5,r K ω j,r + λ 8,r C ω j,r

U ω j,r = U 0,r + j l=1 y ω l,r , j ∈ J, r ∈ R, ∀ω ∈ Ω, (5m)

I ω j,r ≥ 0 j ∈ J, r ∈ R, ∀ω ∈ Ω, (5o)

K ω j,r ≥ λ 6,r H ω j,r − (U ω j,r − C ω j,r ) j ∈ J, r ∈ R, ∀ω ∈ Ω, (5p)

j∈J r∈R y ω j,r e 1 ≤ ∆ ∀ω ∈ Ω, (5r)

y ω t(n),r − y n,r = 0, z ω t(n) − z n = 0, η ω t(n) − η n = 0, ∀ω ∈ β(n), ∀n ∈ N,

S ω j,r , I ω j,r , T ω j,r , H ω j,r , C ω j,r , R ω j,r , F ω j,r , B ω j,r , C ω j,r , O ω j,r ≥ 0, j ∈ J, r ∈ R, ∀ω ∈ Ω,

y ω j,r ∈ {0, 1, 2, . . . ,

Objective Function (5a). The objective function (5a) minimizes the total expected number of tested infected individuals and deaths and the conditional value-at-risk over all stages j and scenarios ω.

Population Infection Dynamics Constraints (5b) -(5j). Constraints (5c) give the initial number of the susceptible, tested infected, untested asymptomatic infected, hospitalized, ICU patients, recovered, and deceased individuals, respectively, in each region r at the beginning of the planning horizon. Constraint (5d) makes sure that the number of susceptible individuals in region r at stage j + 1 under scenario ω equals the number of susceptible individuals at stage j minus the number of susceptible individuals who become either tested infected or untested asymptomatic infected at stage j. In this equation, the number of untested asymptomatic infections equals the number of tested infections multiplied by the proportion of the untested asymptomatic infections to the tested infections. Constraint (5e) shows that the number of tested infected individuals in region r at stage j + 1 under scenario ω equals the number of tested infected individuals at stage j plus the infected individuals caused by short-term migration, plus the newly tested infections at time j, minus the recovered and deceased infections of tested individuals at stage j, minus the hospitalized individuals at stage j. Constraint (5f) implies that the number of untested asymptomatic infections in region r at stage j + 1 under scenario ω equals the number of untested asymptomatic infections at stage j plus new, untested asymptomatic infections, minus the recovered untested asymptomatic infections at stage j. Constraint (5g) shows that the hospitalized individuals in region r at stage j + 1 under scenario ω equals the number of hospitalized individuals at stage j plus the newly hospitalized individuals at stage j, minus the recovered and deceased individuals at stage j, minus the individuals who move to the intensive care unit (ICU) at stage j. Constraint (5h) indicates that the total number of individuals in ICU in region r at stage j + 1 under scenario ω equals the total number of individuals in ICU at stage j plus the individuals who moved to an ICU at stage j, minus the individuals who are recovered or died at the ICU at stage j. Constraint (5i) shows that the number of recovered individuals in region r at stage j under scenario ω equals the number of recovered individuals at stage j plus the recovered individuals from tested infected, untested asymptomatic infected and hospitalized individuals, and ICU patients at stage j. Constraint (5j) indicates that the number of deceased individuals in region r at stage j + 1 under scenario ω equals the number of deceased individuals at stage j plus the deceased individuals from tested infected and hospitalized individuals and ICU patients at stage j.

Ventilator Logistics and Capacity Constraints (5k) -(5r). Both (5k) and (5l) are, respectively, equivalent linear constraints of

C ω j,r = min{λ 6,r H ω j,r , U ω j,r − C ω j,r } j ∈ J, r ∈ R, ∀ω ∈ Ω,

with additional linearization variables, using the method presented in Yin and Büyüktahtakın (2021a) ; Kıbış and Büyüktahtakın (2017) . Constraints (5k) [or non-linear equivalent (6)] ensure that the number of individuals admitted to the hospital in region r at stage j under scenario ω equals the minimum number of individuals who require hospitalization and the available hospital capacity at stage j. Constraints (5l) [or non-linear equivalent (7)] imply that the number of individuals admitted to an ICU in region r at stage j under scenario ω equals the minimum number of individuals who require treatment in ICU and the number of available ventilators at stage j. Constraint (5m) represents that the cumulative number of ICU beds (equivalent to ventilators) in region r at stage j under scenario ω equals the initial number of ICU beds plus the cumulative number of ICU beds (new ventilators) allocated from stage 1 to stage j. Constraints (5n) - (5q) show that the number of individuals who can not be admitted to the hospital or the ICU due to limited capacity should be greater than or equal to zero. Constraint (5r) represents that the cost of ventilators allocated over all regions and time stages under scenario ω cannot exceed the total budget allocated for ventilators. The budget here also represents the maximum total ventilator supply that could be available throughout the planning horizon.

Risk Measure Constraints (5s) and (5t). Constraint (5s) indicates the difference between the objective function value and the value-at-risk for each stage j under each scenario ω. Constraint (5t) ensures that the loss value exceeding the value-at-risk is included in the CVaR calculation, and thus z ω j should be greater than or equal to zero.

Non-Anticipativity, Non-Negativity and Integrality Constraints (5u) -(5w). Constraint (5u) is the non-anticipativity constraint, indicating that the scenarios that share the same path up to stage j should also have the same corresponding decisions. Constraint (5v) indicates that the number of individuals in each compartment in region r at stage j under scenario ω should be greater or equal to zero. Constraint (5w) implies that the number of allocated ventilators should be an integer. We implement this mixed-integer linear programming (MILP) formulation (5a)-(5w) for the rest of the paper.

We present a new region sub-problem and lower and upper bounds to reduce the optimality gap of solving our risk-averse multi-stage stochastic programming problem (5a)-(5w). We are inspired by the study of Büyüktahtakın (2021) , where the general risk-averse stochastic MILP problem was decomposed into scenario sub-problems, while keeping all original constraints the same in the sub-problem. However, our sub-problem is different than that of Büyüktahtakın (2021) in that we decompose the problem with respect to the regions rather than scenarios and then derive bounds based on those sub-problems. Further, we decompose the problem into two: region-r and risk sub-problems. The new scenario sub-problem and bounds are described below.

Definition 4.1. (Region-r Sub-problem) The region-r problem (P r ) is formulated as follows:

Specifically, in P r we minimize the objective function value only under region r while keeping all the variables and the constraints from the original problem (5a)-(5w), P .

Proposition 1. (Region-r Lower Bound) Let (I ω j,r * , F ω j,r * ) be the optimal solution for P corresponding to region r. Then the following inequality holds:

(9)

Proof. P and P r share the same feasible region, but P minimizes the full objective function (5a), while P r only minimizes the objective function for only region r (8a). Thus, the optimal solution of P for region r, (I ω j,r * , F ω j,r * ), is feasible but sub-optimal for P r . 2

In order to get a lower bound on the full objective function (5a), we define the risk sub-problem as follows:

Definition 4.2. (Risk Sub-problem) The risk sub-problem is formulated as follows:

Proposition 2. ZR is superadditive onR ⊆ R, which means that for R 1 , R 2 ⊆ R

Proof. P R 1 , P R 2 , and P R 1 ∪R 2 share the same feasible region with P . Thus, the optimal solution of the P R 1 ∪R 2 is feasible but sub-optimal for each of P R 1 and P R 2 . 2 Proposition 3. (Lower Bound) Let Z * represent the objective function of the original problem (5a)-(5w), P . Then we have:

Proof. The proof follows from the generalization of the superadditivity result in Proposition 2 and the equivalence of Z R and r∈R Z r by Definition 4.1. In other words, we have:

2 20

Proposition 4. (Upper Bound) Letẋ r be the optimal solution of region-r problem, P r , and Z(ẋ r ) be the objective value of original problem (5a)-(5w) whereẋ r is substituted in the original problem objective function (5a). Then we have:

Proof. Becauseẋ r is the optimal solution for P r , it is feasible for problem P (5a)-(5w). As Z * ≤ Z(ẋ r ) holds for each r ∈ R, Z * is bounded above by the minimum Z(ẋ r ) over all r ∈ R.2

In our computational experiments, we demonstrate the effectiveness of the proposed region-r bounds (R) in solving our risk-averse multi-stage stochastic programming problem (5a)-(5w). We use the CPLEX solution and the scenario-ω bounds (S) proposed by Büyüktahtakın (2021) as comparison benchmarks under different budget levels with the following solution approaches:

cpx: Direct solution of the model (5a)-(5w) by CPLEX using its default settings. lb: Region-r lower bound inequalities (9) (lb-R) or scenario-ω lower bound inequalities (lb-S). ub: Upper bound inequalities (13) (ub-R) or scenario-ω upper bound inequalities (ub-S). lb+ub: Inequalities (9) and (13) (lb+ub-R) or scenario-ω lower and upper bound inequalities (lb+ub-S).

In Table 3 , we report results for a combination of columns defined as follows.

The total budget level used. Exp: Solution approach. CutTime: Inequality generation time in CPU seconds, including the solution of all subproblems used to generate bounds. Time: Total CPU seconds required to solve the problem, including inequality generation time (CutTime). MIPGap: Final optimality gap for any given Exp (cpx, lb, ub, or lb+ub) at the end of the time limit. GapImp 1 : Percentage improvement in MIPGap using lb, ub, or lb+ub. For example, for ub, (GapImp 1 = 100 × (1-MIPGap ub /MIPGap cpx ), where MIPGap cpx and MIPGap ub are the final optimality gaps by cpx and ub, respectively. IGap: Percentage integrality gap of the formulation before inequalities are added (InitGap = 100 × (bestobj − relaxlb)/bestobj), where relaxlb and bestobj are objective function values of the initial Linear Programming (LP) relaxation and the best feasible solution by cpx, respectively. RGap: Percentage integrality gap of the formulation after inequalities are added (RootGap = 100 × (bestobj − rootlb)/bestobj), where rootlb is the objective function value of the initial LP relaxation after lb and ub are added. GapImp 2 : Percentage improvement in the integrality gap at the root node (GapImp = 100 × (1-RGap/IGap). R: Method using region-r bounds for 5 region sub-problems. S: Method using scenario-ω bounds proposed by Büyüktahtakın (2021) for 40 scenario subproblems. Table 3 presents the results for CutTime, Time, MIPgap, and GapImp 1 over cpx and lb, ub and lb+ub for both R and S. Since a scenario-ω sub-problem is a weaker relaxation of the original problem (5a)-(5w) compared to a region-r sub-problem, lb-S provides a weaker lower bound for the original problem than lb-R, and thus does not improve the Igap given a limited number of scenario sub-problems used in S. Therefore, we only demonstrate Igap, Rgap, and GapImp 2 for R. Each row of Table 3 presents results for an instance with budget levels $30M , $35M , $40M , $45M , and $50M , while the Average row represents average results for all instances. Overall, instances get easier as the budget level increases, GapImp 1 and GapImp 2 due to R show a decreasing and increasing trend, respectively. Both lb and ub for R and S are generated within 10 and 20 CPU minutes, respectively, on average. Our preliminary results have shown that region-r lower bound inequalities (9) performed better than the lower bound (12); thus we do experiments using the inequalities (9). In Table 3 , ub-R provides the best benefit in reducing the MIPgap at the end of the time limit, while lb-R helps the most in reducing the LP relaxation gap. For example, ub-R improves the MIP gap by 15.0%, while lb-R reduces the initial gap by 27.3% on average. Thus, we recommend adding ub-R as a cut into the formulation to reduce the final MIP gap and using lb-R only in the branch and bound tree to strengthen the LP relaxation of the formulation. Note the instance with $50M is relatively easy, and thus the MIP gap is zero after two hours of the time limit, and the LP relaxation gap of 1.1% is completely closed by lb-R. Compared with R, S takes more time to generate and does not provide a better final MIP gap under lb-S and lb+ub-S over the corresponding R bounds. However, ub-S provides a better MIP gap (MIPgap and GapImp 1 ) than ub-R, except for the easiest instance with the $50M budget level.

This section provides the data used to calibrate model parameters and formulate the model, including population and short-term migration data, transmission parameters, as well as the cost of a ventilator. As shown in Figure 4 , we select eight counties that are most impacted by the pandemic in the states of New York and New Jersey for our case study. They are New York County, Kings County, Queens County, Bronx County, Richmond County, Hudson County, Bergen County, and Essex County. In our multi-stage model, each stage represents a two-week period. Thus, all the data regarding transmission and migration are bi-weekly. Table 4 shows the population data for each considered county in New York and New Jersey. Population data is obtained from JHU (2020). The migration rates among the considered counties, estimated from the data on CENSUS (2020), are presented in Table 5 . The blank areas in Table 5 represent a zero short-term migration because the movement among those counties is too small to make an impact on the model results. Table 6 presents the data values for transmission parameters for the studied counties in New York and New Jersey. The data contains the proportion of untested asymptomatic infections, recovery rate, and the death rate for tested infections, hospitalized infections, and ICU patients. Table 7 shows the transmission rate of each county at the first two stages and the impacts of applying different intervention strategies, as discussed in Section 3.3. Table 8 shows the initial number of infections, hospital capacity, and ICU capacity for each county. The data is obtained from JHU (2020). Ventilator Cost. The cost of each ventilator ranges from $5,000 to $50,000 (Glass, 2020). In our case, we consider a cost of $5,000 for each ventilator, and different budget levels are set to impose different upper bounds on the ventilator supply.

This section presents the validation results of the mathematical model in Equations (5a)-(5w) as presented in Section 3 for the 8-stage time period from April 3, 2020, to July 10, 2020. We consider a medium realization of the uncertain asymptomatic proportions at each stage of the planning horizon and compare the number of infections forecasted by our model to the real outbreak data.

The government applied a lockdown strategy from April 3, 2020, to July 10, 2020, at those considered locations in New York and New Jersey. Thus, we use the lockdown strategy and the corresponding transmission rate at each stage in our model for validation. Figure 5 : Comparison of predicted cases with real outbreak data for new infections in New York and New Jersey Figure 5 shows the comparison between the predicted infections and real outbreak data. The model's predictions give a good approximation for the actual number of new infections in each region, implying that the model can capture the disease transmission dynamics under a lockdown intervention strategy. Similar to the validation approaches of Kıbış et al. (2021) and Kıbış and Büyüktahtakın (2019) , we also perform a paired-t-test to analyze the difference between the pairs of predicted new infections and the actual data in each period. As shown in Table 9 , all p-values are greater than 0.05, and thus our model provides statistically similar predictions with the real outbreak data from April 3, 2020, to July 10, 2020, for each considered county. 

We apply our model described in Section 3 to the selected counties in New York and New Jersey. We first solve the risk-neutral model. We incorporate the uncertainty of the proportion of untested asymptomatic infections as well as the short-term migration in the disease transmission and forecast the number of new infections, deceased individuals, hospitalized individuals, and ICU patients under different intervention strategies. Also, we solve the model to determine the optimal location and number of ventilators allocated under different budget levels and scenarios to provide insights into resource allocation over multiple jurisdictions under uncertainty. Due to the high complexity of the mathematical formulation, we solve it for a 5-stage time period. Each stage corresponds to two weeks, resulting in a planning horizon of ten weeks from March 20, 2020, to May 29, 2020. Because each node of the scenario tree has three possible realizations of the random parameters, we solved 243 (3 5 ) scenarios simultaneously.

The mathematical model is solved using CPLEX 12.7.1 on a desktop computer running with Intel i7 CPU and 64.0 GB of memory. We set the time limit at 7200 CPU seconds to solve each test instance. We extend running time for specific budget levels ($30 Million) and interventions ("Lockdown") due to the large optimality gap. In the following subsections, we present results from solving the multi-stage stochastic epidemics-ventilator-logistics model with an application to the case of COVID-19 using the data presented in Section 5.

Here, we present results of our formulation for each time period under different intervention strategies: No intervention ("None"), mask and social distancing for all stages ("Mask and Social Distance"), lockdown for all stages ("Lockdown"), mask and social distancing for the first three stages and lockdown for the following two stages ("Mask + Lockdown"), lockdown for the first three stages and mask and social distancing for the following two stages ("Lockdown + Mask"). The model is solved under the $30 million budget level. The model with the "Lockdown" strategy gives a 4.54% optimality gap after a run time of 43,241 CPU seconds, while the model solved for all other strategies has a zero optimality gap within 7,200 CPU seconds. Figure 6 presents the number of infections and deceased individuals at each stage under different intervention strategies. According to the results, short-term migration influences the number of new infections even under constant transmission rates. As in the first stage, with the same initial transmission rate, the number of infections under different intervention strategies is different from each other. When the stage increases, the difference in the number of new infections among each intervention strategy becomes more and more significant. The "None" strategy has the most infections at each stage, followed by the "Mask and Social Distance" strategy. The "Lockdown" strategy results in the lowest number of new infections compared to those under other strategies at each stage. The "Lockdown" strategy provides a marginally higher number of infections compared to the actual infection data since our model slightly (but statistically insignificantly) overestimates the number of infections. Compared to the "Mask + Lockdown" strategy, the "Lockdown + Mask" intervention leads to fewer infections. This implies that applying the "Lockdown" strategy immediately at the onset of the pandemic followed by the "Mask and Social Distance" intervention is a better strategy than enforcing "Mask and Social Distance" first and delaying the lockdown.

The intervention strategy does not influence the number of deceased individuals as quickly as it impacts the number of infections, as shown in Figure 6 . Here, the number of deceased individuals at the first two stages is much lower than that of the last three stages. Starting from stage three, the number of deceased individuals under different intervention strategies shows a similar trend with the number of new infections. The influence of interventions is further delayed for those confirmed infections to be treated in the hospital and ICU (Figure 7) .

To reduce both the number of new infections and deaths, "Lockdown" is the best strategy. As shown in Figure 6b , the "Lockdown" strategy with the optimal ventilator allocation further reduces the actual number of deaths. Due to the negative impact of COVID-19 on employment and its economic burden, governments are often forced to stop the lockdown and reopen businesses. In such cases, applying "Mask and Social Distance" after a certain period of "Lockdown" will be the best choice. Figure 7 shows the number of hospitalized individuals and ICU patients at each stage under different intervention strategies. Similar to the number of deceased individuals, there are delays in the impact of government interventions on the number of hospitalized individuals and ICU patients. An infected person may have mild symptoms for about one week, then worsen rapidly (Harvard Medical School (2021)). Thus, it may take some time for patients to be admitted to the ICU, so the impact of interventions on the number of ICU patients is delayed one more stage compared to the hospitalized cases. As shown in Figure 7 , a higher number of hospitalized individuals at stage j will lead to more ICU patients at stage j + 1.

For all the stages, the "Lockdown" strategy has the lowest number of hospitalized individuals and ICU patients, followed by the "Lockdown + Mask" intervention. The ICU patients of "None," "Mask and Social Distance," and "Mask + Lockdown" are almost the same at stages three to five. This is because under those, the need for ventilators is large, and the number of treated individuals in ICUs depends on the minimum number of patients who require ventilators and the ventilator supply in those ICUs. Therefore, the number of treated patients in ICUs is limited by the tight ventilator availability. Table 10 shows the number of ventilators allocated to each region at stages one and two and the total number of ventilators allocated throughout the planning horizon under different budget levels and three select scenarios. The "All Low," "All Medium," and "All High" scenarios represent low, medium, and high realization of the proportion of untested asymptomatic infections at each stage of a five-stage planning horizon, respectively. To analyze the impact of budget on the optimal ventilator allocation decisions, we select $10M as the limited budget level, $20M as the medium budget level, and $30M as the ample budget level. The model has zero optimality gap under the $10M budget level, 0.03% optimality gap under the $20M budget level, and 1.25% optimality gap under the $30M budget level within two hours of solution time.

The results in Table 10 demonstrate that the location and number of ventilators allocated depend on several factors, including the initial and evolving disease transmission rates, the population and the number of initial infections in a region, and the existing ventilator capacity. Thus, the optimal ventilator allocation should be determined case-by-case. more ventilators as the disease transmission rates get higher. Independent from the scenarios, the majority of the regions get more ventilators allocated when the budget level increases from limited to medium and ample. This is in particular observed for the areas with high initial transmission rates (e.g., Kings). The benefit of giving resources to those regions is higher than those regions with low initial disease transmission when the budget increases. Moreover, the model is forced to make difficult decisions, and some of the regions may not have any ventilator allocated under a limited ventilator supply or "All High" transmission scenario under a very limited budget (see, e.g., New York and Hudson counties). When the budget is too tight, the regions either with a high transmission rate or a low capacity get the priority for ventilator allocation. As the budget increases to medium and ample, the model allocates more capacity to the regions with a higher population and a larger initial number of infections but with a lower transmission rate. Also, the stage-wise distribution of ventilators has a high relationship with the available budget. If the budget is tight (see total budgets of $10M and $20M ), most ventilators are distributed within the first two stages. As we increase the budget, some of the ventilators are allocated in stages three and four in addition to stages one and two. Thus, a higher budget level also provides some flexibility in delaying the ventilator allocation to some regions.

In addition to considering a cost of $5, 000 for each ventilator, we perform a sensitivity analysis on different ventilator costs and their impact on the number of ventilators allocated to each region. Specifically, we change the ventilator unit cost from $5, 000 to $10, 000, $15, 000, and $20, 000 under the budget levels of $10M , $20M , and $30M (see Tables S1, S2, and S3 in the Online Supplement S1 for more detailed results). Results show that all ventilators are allocated in time stages one and two when the unit ventilator cost of $5, 000 is increased to $10, 000, $15, 000, and $20, 000. According to the results, the ventilator allocation decisions show similar trends under different ventilator costs. For instance, when the ventilator cost increases and the disease transmission situation worsens (medium and high scenarios), more ventilators are allocated to Bergen County. When the ventilator cost increases from $5, 000 to $10, 000, $15, 000, and $20, 000, New York County and Kings county receive almost the same number of ventilators allocated under different budgets and scenarios. Under the "All High" scenario and $30M budget level, Bronx County receives more ventilators than in other scenarios due to the ample budget.

Furthermore, we can find some trends in the ventilator allocations under different costs. Independent from the budget level and scenario, when the ventilator cost increases, the number of ventilators allocated to Queens and Bronx decreases, while other regions have a similar number of ventilators allocated (e.g., Kings). This is because the model prioritizes allocating an appropriate number of ventilators to these regions with a smaller number of ICUs to support the essential operation of the ICU. When the ventilator cost decreases (or more budget is available), more ventilators are allocated to the regions with high transmission rates and a high initial ICU capacity to reduce deaths.

In this section, we perform an analysis of the risk parameters ϕ and α in terms of their impact on the expected number of infected and deceased people as well as the CVaR of the impact. Specifically, under the $30M budget level, we compare four different problems with respect to their risk-averseness level, adjusting ϕ and α values accordingly-risk-neutral (ϕ = 0, α = 0.95), weak risk-aversion (ϕ = 1, α = 0.3), mild risk-aversion (ϕ = 10, α = 0.6), and strong risk-aversion (ϕ = 10, α = 0.95). The model under the mild risk-aversion results in a high optimality gap (13%) within 7200 CPU seconds running time. Therefore, we solve the scenario-ω problems described in Section 4.1 and obtain the lower and upper bounds for the original problem. For our problem, we select five representative scenarios and add bounds based on the results of those select scenarios in the risk-averse model. After implementing the scenario bounds, the optimality gaps over all of the risk-averseness levels reduce to less than 9.13%.

We decompose the objective function (5a) into the Expected Impact [E(f (x, ω))] and the Expected Risk [CVaR α (f (x, ω) )], as demonstrated in Equation (4), to analyze the impact of risk trade-off on the results. Table 11 presents the value of the objective function (5a), expected impact, and expected risk (without ϕ) under different risk-averseness levels. Specifically, the expected impact represents the expected total number of infections and deceased individuals, and the expected risk corresponds to the expected CVaR term in (5a) without the ϕ value. According to Table 11 , the risk-neutral formulation has the highest expected risk since it does not consider the risks related to high-stake scenarios. When comparing different risk-averseness levels, as both ϕ and α increase, the level of risk-averseness and the expected risk increase. The optimal objective function value increases due to the additional risk term added into the objective formulation. The expected impact also increases, implying the cost of being risk-averse, which is the increased number of infections and deceased individuals while trying to mitigate specific disastrous scenarios. The expected impact and expected risk (without ϕ) for various combinations of ϕ = {0, 1, 10} and α = {0.3, 0.6, 0.95} under the $30M budget level are presented in Table 12 . We observe the change of the expected impact and expected risk when changing one of the risk parameters and fixing all other original values. According to the results, when ϕ = 0, the expected risk has the highest value under each α value (because the risk measure is omitted). Fixing the α value, both expected impact and expected risk show an increasing trend due to the increase of ϕ from 1 to 10. When we move from risk-neutral (ϕ = 0) to risk-averse (ϕ = {1, 10}), the expected impact always increases. Similarly, ϕ = {1, 10} increases the expected impact compared to the risk-neutral model. Besides, when fixing the ϕ value and increasing the α value, the expected risk increases because we increase the confidence level for reducing the risk of having an extremely large number of infections and big losses of lives. In this paper, we present a general multi-stage mean-risk epidemics-ventilator-logistics model and apply this model to control COVID-19 in select counties of New York and New Jersey. We first explicitly formulate the uncertainty of the proportion of untested asymptomatic infections at each stage, generating a multi-stage scenario tree. We then develop a compartmental disease model and integrate the short-term human movement among multiple regions into this model. We also derive a time-and space-varying disease transmission formulation and a logistics sub-model. We then integrate all those components into one mathematical formulation, which minimizes the number of infections and deceased individuals under different intervention strategies. We also present new bounds on the problem and show their computational benefits in terms of the optimality and relaxation gaps. We find that region bounds provide a better improvement on the LP relaxation gap, while scenario bounds provide a better MIP gap for harder instances with the expense of longer sub-problem solution times.

We solve the epidemics-ventilator-logistics model under different budget levels to determine the ventilator-distribution optimal timing and location under various pandemic scenarios. Next, we apply the CVaR in a nested form over a five-stage planning horizon to minimize the total expected number of infections and deceased individuals, as well as the weighted risk of the loss. Finally, we solve the scenario sub-problems under various scenarios to generate the lower and upper bounds for the original problem, reducing the optimality gap. Our results provide key insights into the resourceallocation decisions for controlling COVID-19 and can be adapted to study the transmission and logistics of other similar diseases.

According to the results, the number of infections, deceased individuals, hospitalized individuals, and ICU patients indicates that short-term migration influences the number of infections, even if the transmission rate is constant over time. The impacts of government interventions on the number of deceased individuals, hospitalized individuals, and ICU patients are delayed because deaths and hospitalization have a higher lag period compared to zero or a small lag period in the growth of infections. Furthermore, the number of ICU patients at each time period depends on the minimum number of patients who require the ICU and the available ventilators. Thus, the number of ICU patients might be at the capacity limit even under different intervention strategies at a particular stage. The "Lockdown" strategy is the best way to control disease transmission. Nevertheless, "Mask and Social Distance" applied after the several stages of "Lockdown" is the second-best strategy to alleviate the pandemic's economic impacts.

The region with a high initial transmission rate and low initial ICU capacity receives more ventilators. This is because other regions with low initial transmission rates have fewer infections, even if they have smaller initial ICU capacity. Independent from the transmission scenario, an area with a high initial transmission rate has more ventilators allocated. This is because the benefit of giving resources to those regions is higher than the regions with low initial disease transmission when the budget increases. Moreover, when the budget is limited, all of the ventilators are allocated at the first two stages. When the budget becomes ample, decision-makers would have some flexibility in delaying ventilator allocation to later stages of the planning horizon.

Overall, we can derive general insights based on our sensitivity analysis. However, if the decisionmakers prefer to be precise in the optimal allocation, they are encouraged to solve a complicated optimization model as presented in our paper. This necessity shows the importance of using complex mathematical models to tackle the difficulty of the pandemics control problem over multiple regions and time periods.

The increase in the mean-risk trade-off coefficients in the risk-averse model improves the confidence level, reducing the loss in the right tail of the objective function values (the number of infected and deceased individuals over highly-adverse scenarios). However, we should expect more infections and deceased individuals on average considering all possible scenarios when we want to decrease the impact of adverse scenarios by increasing the risk-averseness level.

This study leads to several future directions for research. For instance, vaccine allocation is also essential as it can potentially protect people from being infected. The combination of vaccine allocation and other interventions will provide more flexible strategies to prevent and control the disease. For example, for the region with a low transmission rate and high vaccine coverage, decision-makers could consider lifting the "Lockdown" earlier to stimulate the economy. Furthermore, the mathematical model cannot allocate ventilators to some regions under a very tight budget, and so future research could investigate ethical and fair resource allocation strategies during a pandemic. Also, some of the assumptions and inferences made in our model could be updated in a future study as more data becomes available. The proposed multi-stage mean-risk epidemics-ventilator-logistics formulation is highly complex and also interesting in terms of developing new solution methodologies. In this paper, we compare our region bounds with the scenario bounds of Büyüktahtakın (2021) . Future research could also investigate the performance of the scenario dominance decomposition cuts of Büyüktahtakın (2021) and other prominent decomposition approaches (e.g., the stochastic dual dynamic integer programming method of Zou et al. (2019) ) proposed for multi-stage stochastic mixed-integer programs.

Risk management for forestry planning under uncertainty in demand and prices

Mathematical modelling to assess the impact of lockdown on COVID-19 transmission in India: Model development and validation

Association between mobility patterns and COVID-19 transmission in the USA: A mathematical modelling study

Sampling scenario set partition dual bounds for multistage stochastic programs

The standard of care of patients with ARDS: Ventilatory settings and rescue therapies for refractory hypoxemia

More COVID-19 patients are surviving ventilators in the ICU

From predictions to prescriptions: A data-driven response to COVID-19

COVID-19 (SARS-CoV-2) ventilator resource management using a network optimization model and predictive system demand

Reallocating and sharing health equipments in sanitary emergency situations: The COVID-19 case in Spain

China's successful control of COVID-19

Risk-averse multi-stage stochastic optimization for surveillance and operations planning of a forest insect infestation

A simulation-deep reinforcement learning (SiRL) approach for epidemic control optimization

Stage-t scenario dominance for risk-averse multi-stage stochastic mixed-integer programs. Accepted for publication in Annals of Operations Research

A new epidemics-logistics model: Insights into controlling the Ebola Virus Disease in West Africa

Demand for hospitalization services for COVID-19 patients in Brazil

CDC, 2021. Clinical Questions about COVID-19: Questions and Answers

Healthcare impact of COVID-19 epidemic in India: A stochastic mathematical model

New Zealand and Australia were Covid success stories. why are they behind on vaccine rollouts?

Stochastic dynamic resource allocation for HIV prevention and treatment: An approximate dynamic programming approach

Epidemics control and logistics operations: A review

Multistage stochastic programming: A scenario tree based approach to planning under uncertainty, in: Decision theory models for applications in artificial intelligence: concepts and solutions

ArcGIS map

The behavioural challenge of the COVID-19 pandemic: Indirect measurements and personalized attitude changing treatments (impact)

A decision support system for demand management in healthcare supply chains considering the epidemic outbreaks: A case study of Coronavirus Disease 2019 (COVID-19)

Symptoms, spread and other essential information about Coronavirus and COVID-19

Watch: Ventilators are in high demand for COVID-19 patients. How do they work?

The risk of COVID-19 transmission in train passengers: An epidemiological and modelling study

Stockpiling ventilators for influenza pandemics

COVID-19 United States Cases by County

Analyzing bioterror response logistics: The case of Smallpox

A time series-based statistical approach for outbreak spread forecasting: Application of COVID-19 in Greece

Optimizing invasive species management: A mixed-integer linear programming approach

Optimizing multi-modal cancer treatment under 3d spatio-temporal tumor growth

A multistage stochastic programming approach to the optimal surveillance and control of the emerald ash borer in cities

Impact of delays on effectiveness of contact tracing strategies for COVID-19: A modelling study

Epidemiological benchmarks of the COVID-19 outbreak control in China after Wuhan's lockdown: A modelling study with an empirical approach

Early dynamics of transmission and control of COVID-19: A mathematical modelling study. The Lancet Infectious Diseases

A flexible method for optimising sharing of healthcare resources and demand in the context of the COVID-19 pandemic

Human mobility trends during the early stage of the COVID-19 pandemic in the United States

Forecasting COVID-19 and analyzing the effect of government interventions

The good, the bad and the ugly of government responses to COVID

COVID-19 transmission trajectories-monitoring the pandemic in the worldwide context

Bounds on risk-averse mixed-integer multi-stage stochastic programming problems with mean-CVaR

A simplified math approach to predict ICU beds and mortality rate for hospital emergency planning under COVID-19 pandemic

The truth about COVID-19 and asymptomatic spread: It's common

A model of supply-chain decisions for resource sharing with an application to ventilator allocation to combat COVID-19

The asymptomatic and pre-symptomatic spread of COVID-19

Risk aversion in multistage stochastic programming: A modeling and algorithmic perspective

Intubation and ventilation amid the COVID-19 outbreak: Wuhan's experience

Forecasting COVID-19 impact on hospital bed-days. ICU-Days, Ventilator-Days and Deaths by US State in the Next

Optimal resource and demand redistribution for healthcare systems under stress from COVID-19

Impacts of epidemic outbreaks on supply chains: Mapping a research agenda amid the COVID-19 pandemic through a structured literature review

Critical supply shortages-the need for ventilators and personal protective equipment during the COVID-19 pandemic

Selectively caring for the most severe COVID-19 patients delays ICU bed shortages more than increasing hospital capacity

Conditional value-at-risk for general loss distributions

A scalable bounding method for multistage stochastic programs

Finding optimal vaccination strategies under parameter uncertainty using stochastic programming

A mathematical model for COVID-19 transmission by using the caputo fractional derivative

Impact of social distancing measures on COVID-19 healthcare demand in Central Texas

Impacts of transportation and meteorological factors on the transmission of COVID-19

Locally informed simulation to predict hospital capacity needs during the COVID-19 pandemic

A framework for rationing ventilators and critical care beds during the COVID-19 pandemic

Covid-19: An agent-based simulation-optimization approach to vaccine center location vaccine allocation problem

A multi-stage stochastic programming approach to epidemic resource allocation with equity considerations. Published online ahead of print in Health Care Management Science

Risk-averse multi-stage stochastic programming to optimizing vaccine allocation and treatment logistics for effective epidemic response

Resource allocation for epidemic control over short time horizons

Mathematical model for Coronavirus Disease 2019 (COVID-19) containing isolation class

Effects of human behaviour changes during the COVID-19 pandemic on influenza spread in Hong Kong

Stochastic dual dynamic integer programming

We gratefully acknowledge the partial support of the National Science Foundation CAREER Award co-funded by the CBET/ENG Environmental Sustainability program and the Division of Mathematical Sciences in MPS/NSF under Grant No. CBET-1554018. We also thank four anonymous referees and the editors, whose remarks helped to improve the content and clarity of our exposition.

According to the results, the total number of ventilators allocated increases in the budget level due to the high need for ventilators. As shown in Table 10 , some regions with many initial infections, e.g., New York County, do not receive ventilators. This situation occurs because the initial ventilator capacity of those regions is higher than in other counties. Also, results suggest that more ventilators should be distributed to other areas with a higher initial transmission rate than New York County, such as Kings and Queens. Both Kings and Queens have higher initial transmission rates and lower initial ICU capacity than New York. Thus, these regions get more ventilators allocated. Also, regions with a smaller population, such as Bergen and Essex in New Jersey, get a large share of ventilators due to their high transmission rates at the beginning of the pandemic and under the "All High" scenario compared to other New Jersey counties, such as Richmond and Hudson.Under a limited budget level, some regions with low initial infections and low initial ICU capacity (e.g., Bronx) get more ventilators allocated when the scenario changes from "All Low" to "All Medium" and "All High." Also, the number of ventilators allocated in regions with a high initial transmission rate (e.g., Kings and Queens) do not significantly increase as the scenario changes from "All Low" to "All High." These regions usually have much more initial ICU capacity for the treatment of their patients because of their large population. On the contrary, the areas with a lower initial transmission rate but less initial ICU capacity may benefit more if they receive