The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a strain of coronavirus which causes the respiratory illness coronavirus disease 2019 . In 2020, on March 11 th , the World Health Organization (WHO) declared the outbreak of SARS-CoV-2 a pandemic [41] . Due to the novelty of the virus, there was (and, at the time of submitting this manuscript, still is) significant uncertainty concerning the severity and mortality of COVID-19. Furthermore, as of October 2020, no vaccine has completed the trials necessary for approving widespread use [24] . Therefore, many countries are enforcing nonpharmaceutical countermeasures [19, 20] , e.g. 1) social distancing, 2) increased public hygiene, 3) travel restrictions, 4) self-isolation (quarantine), and 5) population-wide mass testing for SARS-CoV-2 infection. However, enforcing these countermeasures for long periods of time can have severe economic and social consequences, both at the national and the global scale [28] . Therefore, there is a need for identifying economic strategies for simultaneously relaxing the countermeasures and containing the pandemic.
Model-based decision support systems can be used for exactly this purpose. They use predictive models to assist decision makers in identifying and evaluating candidate strategies (e.g. [1] ). In particular, given a dynamical model of the spread of SARS-CoV-2, economically optimal (open-loop) mitigation strategies can be identified by solving optimal control problems (OCPs) over several months or even years. A key advantage of this approach is that it can directly account for constraints, e.g. related to the capacity of public healthcare systems. However, given the uncertainty surrounding SARS-CoV-2 and COVID-19, it is advisable to implement the optimal mitigation strategies in closed-loop, i.e. to repeatedly update the strategies when new data becomes available. This is referred to as model predictive control (MPC) [32] and is an established method for advanced process control [31] . The predictive capabilities of the underlying model are crucial for the efficacy of the resulting mitigation strategy, and a common challenge is to identify suitable model parameters.
Epidemics are often modelled using deterministic compartmental models [18] , e.g. the classical SIR model, where individuals are either susceptible, infectious, or removed, or the SEIR model which, additionally, takes the incubation time into account. Optimal control of compartmental models was already an active research topic before the SARS-CoV-2 pandemic (see [35] for a review). In particular, optimal control of SIR models has been considered, e.g. for arbitrary social interaction models [2] and to identify time-optimal mitigation strategies [3, 16] . Optimal control of more complex models has also been considered. For instance, Fischer et al. [12] consider optimal control of a model with two species, Bussell et al. [6] demonstrate the importance of closed-loop mitigation strategies (i.e. of incorporating feedback), and Watkins et al. [39] consider MPC of stochastic compartmental models. In [11] the authors determine control strategies to maintain hard infection caps in a disease-vector model based on the theory of barriers. This approach, however, exploits the low dimensionality of the model. Application of these techniques to complex compartmental models, therefore, requires model order reduction.
In response to the SARS-CoV-2 pandemic, many researchers have presented optimal control strategies, for instance based on Pontryagin's maximum principle (e.g. [21, 30, 42] ). These strategies typically involve 1) extended SIR or SEIR models, 2) nonpharmaceutical countermeasures (often social distancing), and 3) minimization of the number of infected as well as the economic cost of the countermeasures (and often other quantities as well, e.g. the number of deaths). Furthermore, they rarely satisfy hard constraints, for instance related to health care or testing capacities. In the following, we highlight some of the key developments in decision support for SARS-CoV-2 mitigation based on optimal control. Gondim and Machado [14] use a model with three age groups to compute optimal quarantine strategies (for susceptible individuals) which minimize the number of infected and the cost of quarantining. Bonnans and Gianatti [4] compute social distancing strategies based on a model with a continuous age structure. Here, the strategies minimize a combination of 1) the number of deaths, 2) the peak number of hospitalized, and 3) the cost of social distancing. Similarly, Richard et al. [33] present optimal social distancing strategies based on a model with a continuous age and infection duration structure, which minimize the number of deaths and the cost of social distancing. Morato et al. [25] compute on-off (also called bang-bang) social distancing strategies which minimize 1) the number of symptomatic infectious people and 2) the duration of the social distancing policies, subject to constraints on intensive care unit (ICU) occupancy. They use extended SIR models. Carli et al. [7] use MPC to compute social distancing and travel restriction strategies for an extended multi-region SIR model, minimizing the cost of the countermeasures and preventing an overload on the hospitals. Köhler et al. [22] use MPC to minimize the number of fatalities caused by COVID-19, subject to constraints on the economic cost of social distancing. They take a modified SIDARTHE model [13] as basis and use interval arithmetic in the MPC to propagate model uncertainties. Finally, Tsay et al. [37] use MPC to minimize the cost of social distancing and testing, subject to an upper bound on the peak number of infectious people who have been tested positive. They use the unscented Kalman filter to estimate the noisy state variables of an extended SEIR model.
In this work, we address some of the key questions that decision makers involved in the mitigation of the SARS-CoV-2 pandemic are facing: 1) Is mass testing alone sufficient to avoid overloading of ICUs? 2) If not, how much social distancing is then required? 3) How much can social distancing measures be reduced by targeting specific age groups? 4) How do strategies obtained by short and long-term planning differ? 5) What are the benefits of increasing the daily testing capacity or the ICU capacity? Here, the limited ICU capacity is considered as an example for constraints imposed by the health care system or political considerations. Of course, different constraints such as limited personnel for contact tracing could be incorporated as well.
We address the above questions by proposing a novel compartmental model and using optimal control as well as MPC to compute open and closed-loop social distancing and testing strategies. The model contains three age groups, and it accounts for several of the key challenging characteristics of COVID-19, i.e. 1) the incubation time, 2) different levels of symptom severity depending on age, 3) delay of testing results (and the following self-isolation), and 4) delay of hospitalization. Furthermore, we choose values of the epidemiological model parameters based on the current state of knowledge in order to ensure that our numerical results match reality. For concreteness, we use the COVID-19 outbreak in Germany to determine parameters depending on demographics and the health care system. However, we expect our conclusions to carry over to outbreaks in other developed countries as well.
The remainder of this paper is structured as follows. In Section 2, we describe the novel compartmental model of the SARS-CoV-2 outbreak in Germany, and in Section 3, we motivate our choice of model parameters. In Section 4, we demonstrate that optimal control can be used as a decision support tool based on the proposed model, and we conclude the paper in Section 5.
In this section, we propose a dynamical model tailored to COVID-19. The aim is to be able to evaluate the effect of population-wide mass testing (in combination with quarantine) and social distancing measures on the development of the pandemic. To this end, we extend the well-known SIR model.
We start with an illustration of the connection between 1) infectious disease models based on randomly acting individual agents and 2) their approximation by ordinary differential equation compartmental models. This exposition will highlight the interpretation and conversion of parameters when moving from a random to a deterministic model. For simplicity, we consider the classical SIR model in this subsection. However, the connection, especially the interpretation of parameters, is similar for more complex models such as the one described in Section 2.2.
Consider a population of n pop individuals or agents each being either susceptible, infectious or removed. At time t ∈ [0, ∞) denote the (random) set of susceptibles by S t , the set of infectious by I t and the set of removed by R t . Time is modeled continuously and measured in days.
We assume a homogeneous population with contacts between agents a and b following a Poisson process with intensity λ which does not depend on the agents considered. Infections occur randomly upon contact with a fixed probability α if one of the agents is susceptible and the other infectious. Thus, potentially infectious contacts also follow a Poisson process with respective intensity αλ. Similarly, we model other events, in this simple model only recoveries, to occur according to a Poisson process. This implies that the time an agent spends in the infectious compartment is exponentially distributed with rate η, say, which we also assume to be the same for each agent (see [27] for models where these quantities follow other distributions).
We denote by S(t) = E |St| npop , I(t) = E |It| npop and R(t) = E |Rt| npop the expected share of the population which are susceptible, infectious and removed, respectively. Since for large n pop the change of |St| npop over a short time interval can, due to the law of large numbers, be well approximated by its expectation, S(t) will provide a sufficient approximation of |St| npop over the finite time horizon considered for a country the size of Germany. By the same argument, I(t) and R(t) approximate |It| npop and |Rt| npop , respectively, sufficiently well. If a is susceptible he will transition to the infectious compartment upon having an infectious contact. At a fixed time t with a ∈ S t , there are two possible sources of infection for a: either some b ∈ I t which is already infectious or some c ∈ S t which will become infectious himself at some later time.
To determine the probability that b infects a in the time frame (t, t + ∆t], we analyze two competing events: The first is an infectious contact between a and b, and the second is b's recovery from the infectious state. Both events happen independently of one another with exponentially distributed time of occurrence, the first with rate αλ and the second with rate η. Thus the first time of occurrence of one of these is again exponentially distributed with rate αλ+η and the probability that the first occurrence is an infectious contact is αλ αλ+η . In total
For c ∈ S t to infect a in (t, t + ∆t], c has to become infectious himself before he in turn can infect a. This happens only with probability o((∆t) 2 ) and can, thus, be neglected in the following calculations. In total a is moved out of the susceptible compartment with probability
Approximating |St| npop and |It| npop by S(t) and I(t) using the law of total expectation yields
As we assume the time from infection to removal to be exponentially distributed with rate η, a similar but more straight-forward calculation reveals
where η −1 is the mean stay of a single agent in the infectious compartment. We now set β = n pop αλ, which can be interpreted in this model as the daily amount of (potentially) infectious contacts a single agent has. Since S(t) + I(t) + R(t) = 1 for all t, we obtain the following system of ODEs:
To determine suitable parameter values for β and η in this simplistic model, we reiterate that these are best thought of in the probabilistic setting. For the coefficients of the linear terms on the right-hand side, the interpretation is straight-forward: it is the rate of the exponential distribution underlying the time until an agent leaves the compartment. Its inverse is the mean stay in this compartment. For coefficients of interaction (product) terms, here β, the interpretation is the rate at which an agent in the first compartment causes other agents to leave the second compartment. In our setting, this is the daily amount of infections one infectious agent causes which can readily be seen from the definition of β. See Section 3 for a more detailed discussion of the parameter values we use in our model. As above mentioned, these interpretations for the parameters carry over in a straight-forward manner to more sophisticated models such as the one considered in the following.
The SIR model provides a good starting point to study the dynamics of pandemics. However, due to its simple structure it is not suited to model the COVID-19 pandemic adequately. In particular it does not include hospitalization, age-specific disease progressions and interventions. Therefore, we extend the SIR model in three ways.
1. We introduce eight additional compartments. In detail, we take into account that people can be infected, but not yet infectious. We call them exposed (or latent) and denote the compartment by E, see also [18] . Moreover, we split the infectious compartment into three depending on how the course of the infection will be. We distinguish between severe cases I S that are going to need intensive care, i.e. they will move to H ICU at some point in time; mild cases I M that are going to visit a physician and hence be quarantined, i.e. removed; and asymptomatic cases I A that might recover without being detected. Furthermore, we incorporate the possibility of being tested but not yet detected by introducing the compartments T S (severe) and T O (other). We assume that the patients with severe cases will visit a physician at some point before being sent to an ICU. To this end, we introduce P as a pre-ICU compartment which comprises isolated patients at home or on a regular hospital ward. Moreover, we split the compartment of removed people into known and unknown cases R = R U + R K .
2. Each compartment is further divided into n g groups, n g ∈ N, depending on the age of a subject in order to study how these groups affect each other.
3. Social distancing and hygiene measures affect the contact rate as well as the transmission probability. Therefore, β can be used as a time-dependent control input β(t).
The resulting SEITPHR model reads aṡ
where the subscript i ∈ {1, 2, . . . , n g } denotes the age group in ascending order. We enforce ng i=1 N i = 1, where N i denotes the relative size of age group i. We assume a mean incubation time γ −1 independent of both the course of infection and the age of the patient. However, depending on the age, there are different probabilities π S i , π M i , and π A i for the three courses of infection, where π S i + π M i + π A i = 1 for all i. Similar to the SIR model (1), the parameters denoted by η correspond to people being removed from the system, i.e. η S and η M denote those who visit a physician and, therefore, are put into quarantine immediately, while η A represents unreported recovery. We denote the total number of susceptibles and unreported cases by
describes the rate of those being tested per day, where tests are distributed uniformly at random among all individuals in U i . In addition, symptomatic cases who visit physicians are assumed to be tested as well. Therefore, the total number of tests at time t ≥ 0 is given by
Note that testing does not affect the state of non-infectious subjects. Parameters τ S and τ O denote the rate from being tested positive to being detected, and hence being put into quarantine. Furthermore, ρ is the rate from pre-hospital quarantine to hospitalization and σ from hospitalization to being reportedly removed, i.e. σ incorporates both mortality and recovery rate of hospitalized patients. The basic structure of the SEITPHR model (2) is depicted in Figure 1 . Figure 1 : Flow of the SEITPHR model for one age group. The controls are indicated with dashed red edges. Unreported compartments are highlighted by the left red triangle, while tested and detected compartments are highlighted by the right blue trapezium.
For a concise notation we stack the state vectors into x = (x 1 , . . . x ng ) and the controls into u = (β, θ), where
and β = (β ij ) ng i,j=1 with β ij : R ≥0 → R ≥0 and θ = (θ 1 , . . . , θ ng ). Similarly, we denote π ∈ R 3ng , τ ∈ R 2 , and η ∈ R 3 . Thus, we write system (2) aṡ
where p = (π, η, τ, ρ, σ, γ) ∈ R 3ng+8 collects all parameters. Furthermore, we introduce the initial condition
Before we present our choices for the parameters of model (2), let us reiterate that some of the parameters of our model depend on age. We indicate this dependence by an appropriate index which we drop if the parameter is constant across agegroups. For example, π S i is the age-dependent probability of having a severe course of disease while we assume η A , the rate with which asymptomatic cases recover, to be age-independent. N i We use data on the population size of Germany at the end of 2019 from the GENESIS-Online Database of the DESTATIS [36] . The first age group consists of individuals aged younger than 15 years, the second of those older than 15 but younger than 60 years, and the last comprises all individuals older than 60 years. These groupings result in proportions N 1 = 0.14, N 2 = 0.58 and N 3 = 0.28.
The rate at which an infected agent in compartment I j infects susceptibles in compartment S i depends on the contact structure of a population as well as the probability that a contact between a susceptible and infectious agent leads to a transmission of the disease. We base our contact process on data from the POLYMOD study on daily contacts in several European countries [26] . From this data we calculate a contact matrix C = (c ij ) whose (i, j)-th entry is the mean amount of contacts an individual in age group i has with age group j; here we only consider those contacts labeled as physical, since those are more likely to lead to viral transmission.
Let us denote by β 0 ij the rate at which a single infectious agent from age group j infects susceptible agents from age group i if no countermeasures, such as social distancing, are in place. We model β 0 ij to be proportional to c ij and let α be the corresponding proportionality constant. If a single infectious agent is introduced without interventions such as test, quarantine and social distancing measures in place into the otherwise completely susceptible, i.e. virgin, population, the mean amount of secondary cases he causes is the basic reproduction number
There is a wide variety of estimates for R 0 in the literature [29] , with most estimates in the interval [2, 3.5] . We choose a value of R 0 = 2.5 as early, higher estimates might be biased upwards due to imported and undetected cases. Fixing η −1 = 6 (see the discussion on η A below) we calculate α = 5.79% and in turn β 0 ij from (4):
γ The rate at which latent cases become infectious is the inverse of the mean incubation time. This parameter is modeled age-independent and chosen to be 0.19, which corresponds to a mean incubation time of 5.2 days [23] .
These parameters denote the proportion of individuals in age group i that have severe, mild or asymptomatic course of disease. For Germany, the Robert Koch Institut (RKI) has published data on severity of disease progression for 12, 178 cases by age-groups [34] . For our purposes we define a severe case to be a case that will eventually be admitted to intensive care, a mild case being one developing influenza-like symptoms, pneumonia or being admitted to hospital for other reasons. All other cases we classified as asymptomatic.
Thus we obtain π i = (π S i , π M i , π A i ), the proportion of severe, mild and asymptomatic cases in age group i, respectively, as (π 1 , π 2 , π 3 ) = 1 100 Observe that the oldest age group is at highest risk with 3.02% of infected individuals admitted to ICU. Also the proportion of severe cases in the youngest age group is higher than in the middle age group. This might be explained by the fact that cases in the youngest age group are detected less frequently due to them being tested less, leading to overreporting of severe cases. η S , η M , η A These are the rates at which infectious individuals are removed from the infection process, if no mass-testing is implemented, i.e. if θ i = 0. For individuals with severe or mild course of disease when they develop symptoms leading to self-isolation, quarantine prescribed by a physician, or to direct hospitalization. One characteristic of COVID-19 is that even presymptomatic cases transmit SARS-CoV2 [17] . We assume the time from being infectious to symptom onset to be 2 days after which we add 2 more days which it takes before the infectee visits a physician. Thus we choose η S = 0.25.
For mild progressions we assume the same mean duration from being infectious to symptomatic, though in this case individuals self-isolate, visit a physician or receive a positive test result after a mean waiting time of again 2 more days; consequently, we also set η M = 0.25.
For asymptomatic cases in I A i the only way to be removed from the infection process is by recovery from the infection. In [40] positive virus samples were found in patient's throats for up to 8 days after symptom onset. Assuming a lower viral load for asymptomatic cases with only 4 days of potential infectiousness and adding the 2 days of presymptomatic transmission we chose η A = 0.17, corresponding to a mean time of 6 days to recovery for asymptomatic cases. τ O , τ S As we assume the testing related to the controls θ i to be of a random nature, tested individuals are not yet removed from the infection process. Instead we assume positive test results to become available after a mean delay of 2 days. However, severe cases may visit a physician and thus go into immediate quarantine before receiving their test result. The latter transition occurs with rate η S , and hence the faster transition occurs with rate τ S = 1 2 + η S = 0.75. Non-severe cases that are tested, T O i , are removed if they recover naturally (with rate η A ), or receive a positive test result, or visit a physician. That leads to τ O = η M + η A + 1 2 = 0.92 for each age group.
ρ This parameter is the rate at which severe cases move from being in the pre-ICU state to the intensive care unit. This includes time spent in quarantine at home as well as time spent in the hospital in normal care while being isolated.
In [10] the median time from symptom onset to being in intensive care for 50 patients was 9 days. As the median of an exponential distribution Exp(ζ) is log 2 ζ we choose a mean stay of 9 log 2 − 2 ≈ 10.98 = ρ −1 days, accounting for the mean two days from symptom onset to the transition into the pre-ICU compartment.
σ This is the mean time spent on intensive care until discharge or death. According to [10] patients with acute respiratory distress symptoms (ARDS) spent a median amount of 13 days in intensive care and patients without ARDS spent a median amount of 2 days in intensive care. Of the 50 patients considered in this study, 24 were afflicted with ARDS. Converting again between median and mean for the assumed exponential distribution yields a mean time of σ −1 = The DIVI-Intensivregister offers daily information on the amount of free intensive care hospital beds in Germany. On 20 October 2020 they reported a capacity of 8, 872 free beds with 879 actively treated COVID-19 patients [9] . We therefore round the maximal ICU-beds available for COVID-19 patients to 10, 000.
T max In late August until the beginning of October the RKI conducted between 1 and 1.2 million weekly SARS-CoV-2 tests in Germany. This motivates our upper bound T max = 1,200,000 7 of daily tests.
x 0 We initialize our model at time t = 0 with entries of x 0 set to 0 except for those related to the susceptible, latent and infectious compartments. Our choice of initial values is informed by the number of active cases reported by the RKI in late march assuming the proportion of underreporting to be 50%. We hence set the total number of infectious agents at t = 0 to 524 and the number of latent agents to 1672 distributed among the age-groups according to N i . As we explain in Remark 1 our model is robust against misspecification of the initial values. Figure 2 demonstrates the simulation capabilities of our model. Here, the course of the pandemic is visualized if no countermeasures are implemented, i.e. no social distancing (β(t) = β 0 ) and no mass-testing (θ(t) = 0). As expected, the pandemic evolves too fast to satisfy any reasonable cap on the number of required ICUs. In particular, the number H ICU of required ICUs exceeds 100,000 whereas we noted above that in Germany only about 10,000 ICUs are available to treat COVID-19 patients. Therefore, countermeasures are indispensable to avoid an overload on the hospitals. Note that if we distinguish different age groups the pandemic evolves faster, but less ICUs are required, as the pandemic spreads mostly in the less vulnerable, younger age groups. Similar observations, viz. herd immunity being achieved faster in heterogeneous populations in comparison to homogeneous ones, have already been made by [5] .
In this section, we provide information on how to keep the epidemic manageable.
To this end, we formulate suitable optimal control problems (OCPs) and solve them numerically. Since we neither take vaccines nor re-infections into account, we consider the epidemic to be over once herd immunity is achieved, i.e. a state where the introduction of new infectious agents does not lead to an outbreak. Therefore, our main goal is to reach herd immunity with as few social distancing as possible while maintaining strict limits on the ICU occupancy to avoid a breakdown of the health care system. We call a control u = (β, θ) of the system (3) 
, . . . , n g }, and
is satisfied for all t.
A natural stopping point for simulations is when the share of susceptibles has decreased enough to ensure herd immunity even when all countermeasures are lifted completely. The time-dependent effective reproduction number R(t), that is the mean number of secondary cases a primary cases will cause at time t, can be used to determine whether herd-immunity has been reached: this will be the case if R(t) is less than 1. If there is only one infected compartment, as in a simple SIR model, the latter condition is equivalent toİ(t) < 0. If there is more than one infected compartment, as in model (2), [8, 38] have suggested to compute R NGM (t), based on the so-called next-generation matrix, as a proxy for the effective reproduction number. Then, R NGM (t) exhibits the same threshold property as R(t), that is R NGM (t) < 1 implies herd immunity. Thus we use R NGM (t) to check whether our simulations have reached herd immunity. A time horizon of two years (104 weeks) sufficed for all our simulations.
This section is structured as follows. First, we verify the existence of a feasible testing strategy, i.e. without enforcing social distancing. Note that due to delays in testing, the existence of a solution is not trivial and depends on the initial value. Next, we establish an upper bound on the maximal number of tests per day and investigate to what extent social distancing is required in order to ensure feasibility. Throughout our simulations, we assume the length of one control interval to be one week. This reflects the practical constraint that the government cannot change policies arbitrarily often but more realistically on a weekly basis. Throughout our simulations we use the Matlab-inherent sequential quadratic programming (SQP) tool to solve the OCPs.
Here, our goal is to maintain a hard cap on the number of required ICUs with as few tests as possible without enforcing social distancing, i.e. β ≡ β 0 . To this end, we solve the OCP
The objective function penalizes the total number of tests over the entire time horizon [t 0 , t f ] with t f > t 0 ≥ 0. The equality constraint (6c) captures the system dynamics while the one-sided box constraints (6e) ensure that the testing rates cannot be negative. Figure 3 depicts the optimal controls as well as the total number of tests and the number of detected cases per day while Figure 4 shows the impact on the evolution of the epidemic. Here, we computed the effective reproduction number R NGM (t) at each time step, demonstrating that we reached herd immunity.
We observe that there exists a testing strategy that ensures feasibility, which was not obvious from the outset because of the assumed delays. In particular, the bound (6b) is active once it's reached, i.e. H ICU ≡ H ICU max , and becomes inactive when the number of susceptible people falls below a certain threshold and R NGM (t) < 1 indicating the onset of herd immunity. (6) satisfies the socalled turnpike property [15] . Typically, turnpikes indicate the optimal operating state of a system. These are steady states at which the running costs are minimized. In our example, since we do not penalize the number of required ICUs, the best strategy is to stay at the upper bound while saving tests. Once the objective function value is zero the system leaves the state eventually. In particular, regardless of the initial value, the system is steered towards this optimal operating point. As a consequence, a rough estimation suffices as initial guess for our simulations. A rigorous analysis of these turnpikes, however, is left for future research.
However, these results are only of theoretical interest, since this optimal testing strategy would be prohibitively expensive and might not even be implementable at all. For instance, regarding Figure 3 , one observes that the mean testing rate reaches about 0.5, which corresponds to being tested every two days on average. Moreover, the total number of tests per day required for this approach is more than 12,000,000 on average (over 65 weeks), compared to T max ≈ 170, 000 daily tests which are currently conducted in Germany. Note that, even with this enormous testing effort, the number of detected cases, T S + T O , is rather small since the number of infectious individuals is small compared to the total population.
In conclusion, mass-testing alone currently does not suffice to maintain hospitalization caps in reality. These arguments support the government's decision to introduce additional measures like social distancing and hygiene concepts. However, cheap rapid test kits might change the situation favorably as they could be made widely available, self-administered while giving immediate test results.
In the following subsections, we enforce T max as an upper bound on the amount of daily tests. Under this additional constraint we then determine the minimal amount of social distancing required to reach herd immunity. The success of such measures depends on the acceptance and thus compliance by the general population.
In a first step, we determine an optimal social distancing strategy by penalizing the deviation of β from β 0 equally over all age groups. This might increase acceptance in the general population due to the (perceived) fairness of such measures: everyone is treated equally and contacts are reduced by the same proportion for everyone. In reality such strategies may be hard to conceive as different measures affect the age groups differently, i.e. closing schools and nurserys affects those in the lowest age group the most while leaving the oldest age group unaffected. Nevertheless a mixture of many different non-pharmaceutical measures may be able to achieve such a reduction in contacts.
We introduce a time-varying factor δ = δ(t) describing the amount of social distancing that is implemented. Moreover, we choose to penalize the 2 deviation of this control input from δ = 1 in the objective function in order to smooth the optimal control. For instance, penalizing the 1 deviation yields bang-bang controls, i.e. the optimal solution jumps back and forth between the two extremal options: no contact restrictions and lock down (simulations not shown). Therefore, we determine an optimal homogeneous social-distancing policy by solving min θ,δ
Note that we allow to distribute the tests among the age groups by not fixing θ i , but enforcing (7e) and (7g). Furthermore, we introduce a regularization term with weight κ = 10 −5 . The choice of κ is based on simulations. In contrast to (6), we always find a feasible solution of (7) if the epidemic has not yet evolved too far. More precisely, by choosing δ = 0, which corresponds to a complete lockdown, we are (theoretically) able to stop the spread. Therefore, if the initial number of people with a severe course of infection is sufficiently low, the upper bound on the number of ICUs will not be violated. A highly fluctuating social distancing strategy may lead to low acceptance in the general population, because people have to adapt to new rules every few weeks. Thus before we solve (7) let us have a look at what happens if we consider a constant value for δ over time, i.e. a social distancing strategy without fluctuations. Figure 5 (left) shows that fewer contacts result in a longer time for the epidemic to abate on the one hand, but a lower number of total infections within the considered time horizon on the other hand. Moreover, Figure 5 (middle) visualizes that quite strict social distancing is needed in order to meet the ICU capacities. The maximal value of δ to stay feasible is 0.487, i.e. contacts needed to be more than halved over three years. Furthermore, once we lift the restrictions, see Figure 6 , there might be another outbreak. In particular, the stronger the restrictions were in the beginning, the stronger the second outbreak will be. Therefore, it is essential to establish herd immunity before lifting all restrictions, and to adapt the policy over time.
A visualization of the optimal solution of (7) can be found in Figures 7 and 8 . As mentioned above, the bound on H ICU is not violated. Since the weight κ is chosen sufficiently small, the upper bound on the total number of tests per day is active as long as the upper bound on δ is not. However, note that not all age groups are tested equally. More precisely, only the middle-aged group is tested at all. The reason is that this group is the largest (N 2 > N 1 + N 3 ) and has the highest contact rates (c.f. (5) ) and therefore, contributes more to the spread of the epidemic than the other groups. Furthermore, we observe that the social distancing policy has to be quite strict in the beginning. In particular, min t δ(t) ≈ 0.3 which corresponds to a reduction of average contacts per person by 70%. However, this can be qualitatively compared to the measures taken in Germany starting in mid March 2020 when contacts were reduced by school and restaurant closures as well as other contact restrictions.
In conclusion, social distancing is an effective tool to keep the epidemic manageable. Comparing the results of (7) to the simulations with constant δ we see that a (partial) lockdown appears inevitable. However, our simulations suggest to let the epidemic evolve for a few weeks, then enforce a contact reduction down to approximately 30% for 2-4 weeks before slowly lifting the restrictions over the next 12 months until herd immunity is achieved.
The constraint that contacts are reduced by the same proportion for each age group is restrictive and it is plausible that more efficient solutions exist when contact reductions are distributed differently across age groups. One reason to consider such a strategy is that it may be more efficient at stopping the spread of the epidemic; as mentioned above the middle-age group is the driver of the epidemic while the oldest age group consists of the most vulnerable individuals. In any case, such an age-differentiated social distancing strategy needs to be accepted by the whole population to be successful.
Hence, we improve the social distancing policy computed above by allowing it to depend on age. Given the solution (θ , δ ) of (7), we solve the OCP min θ,β
Here, we use δ to define β min ij = min t δ (t)β nom ij , i.e. the lower bound on β in (8) is the worst case of (7). Therefore, no one is treated worse than when applying homogeneous social distancing. Note that (θ , β) with β = δ β 0 is feasible for problem (8) . As in (7) we penalize testing as soon as β(t) = β 0 holds.
Results for (8) can be found in Figures 9 and 10 , wherē
describes the average number of contacts per person and day in a heterogeneous population. Here, we used κ = 10 −5 . The corresponding value for β 0 isβ 0 = Figure 9 : Optimal age-dependent social-distancing strategy for three age groups over two years. Figure 10 : Evolution of the compartments associated with controls depicted in Figure 9 .
0.4167. This allows to compare the solution β ij (t) withβ 0 δ (t) obtained from (7) . Similar to the solution of (7), the upper bound on testing is active most of the time, while essentially only the middle-aged group is tested. The social distancing measures are less restrictive than for (7) which makes compliance with the measures more likely. However, the measures could be perceived as unfair, since the contacts of the oldest age group are restricted most. Moreover, the contacts of the middle-age group are least restricted. Therefore, the working class would be allowed to go to work which benefits the economy.
In conclusion, social distancing is crucial to avoid an overload on the hospitals. In addition, testing middle-aged people helps to reduce the required amount of social distancing. Furthermore, all presented strategies support a lock down a few weeks into the epidemic, which is followed by lifting the restrictions step by step until herd immunity sets in. Age-differentiated social distancing might be hard to argue for, but it helps to end the epidemic several months earlier and, therefore, support the economy.
The control strategies derived in the previous subsections provide rough guidelines for how the epidemic can be controlled. However, from a decision maker's perspective, it will be hard to argue for policies taking effect in the far future. In particular, there are many uncertainties that might affect the performance of the control strategy over the time span of two years, and hence the control strategy needs to be adjusted over time. Therefore decision need to be revised constantly adapting to the changing conditions during the epidemic's progress. Model predictive control (MPC) provides a state of the art methodology to tackle such Figure 11 : Optimal control for solving (8) in closed loop for varying prediction horizon length. For the sake of readability, we depicted average values of θ and sums of T tot over the age groups. Figure 12 : Evolution of the epidemic based on the controls depicted in Figure 11 . For the sake of readability, only the sum over the age groups is visualized.
problems. The basic idea of MPC is to consecutively solve a series of OCPs over a smaller horizon of K control intervals rather than solving a single OCP over the whole horizon. More precisely, only the first part of the optimal control derived by solving such an auxiliary OCP is implemented. Then, the time window is shifted, and the procedure is repeated based on updated measurements. For a detailed introduction to MPC we refer to [32] . Here, we tackle (8) via MPC; the earlier problems can be treated analogously. The MPC scheme for (8) is summarized in Algorithm 1.
Input: Prediction horizon length K, length of control interval ∆t. Set time t = t 0 . Repeat:
1. Obtain current statesx = x(t).
Results based on varying prediction horizon lengths can be found in Figures 11  and 12 .
The basic structure of both the optimal control and the associated states is comparable to the open-loop solution presented in the previous subsection. Therefore, we stopped the simulations after one and a half years. The length of the prediction horizon affects mainly the optimal social distancing policy. In particular, the larger the prediction horizon, the less social distancing is needed in total. More precisely, for bigger K, we implement a slightly stricter lockdown but can start it later and relax it earlier. Furthermore, the larger K gets the closer the optimal solution is to the open-loop solution. In particular, the MPC solutions qualitatively resemble the open-loop solution: after an early lockdown, social distancing is slowly lifted.
For K = 3, the ICU capacity reaches its upper limit earlier due to the laissezfaire policy in the beginning. However, this constraint also becomes inactive earlier.
For even shorter prediction horizons recursive feasibility cannot be guaranteed, i.e. the ICU cap might be violated (simulations not shown).
So far, we assumed both the upper bounds on the number of tests per day, and on the number of ICUs to be fixed at our chosen values. In practice, these conditions may change: free ICU capacity might exhibit seasonal patterns and the number of possible tests per day depends on infrastructure and available personnel. In addition, varying the upper bounds is useful to illustrate the benefits of increased testing and higher ICU capacities. In this subsection, we investigate the impact of these parameters on the optimal social distancing policy numerically.
First, we study the effect T max has on the social distancing by solving (7) via MPC, see Figure 13 (left). As pointed out in the previous subsection, the prediction horizon length affects the start and end time of measures as well as its peak (simulations not shown). In addition, increasing T max by some factor T max fac ≥ 1 shifts the whole δ curve upwards, i.e., as expected, the more tests are available, the less social distancing is required. Furthermore, Figure 13 (left) visualizes the impact of T max on the objective function value of (7). Figure 13 : Impact of T max on social-distancing costs (left) and of H max on both social-distancing costs (middle) and testing (right). In the last two subfigures the currently available number of ICUs in Germany is highlighted by a vertical dashed line. The dashed horizontal lines in the right-most figure indicate the total testing capacities over the entire simulation horizon. Factor of modification of T max denoted by T max fac .
Second, we investigate the impact of the number of ICUs on the optimal solution of (8) . Results can be found in Figures 13 and 14 . For the simulations in Figure 13 (middle and right) we used MPC with prediction horizon K = 12 weeks. Figure 13 (middle) clearly shows that the number of available ICUs directly affects the cost function value. While for a small value of H ICU max , every additional ICU contributes, for large values, a saturation seems to take place. In particular, doubling the current number of available ICUs does help, but the benefit becomes negligible when increasing it further. These phenomena are almost unaffected by doubling or halving T max . However, when there are not enough ICUs, then the upper bound on T tot is always active, see Figure 13 (right), where T tot is at its maximum value all the time. Moreover, an increase in the number of ICUs clearly Figure 14 : Impact of the available number H ICU max of ICUs and the prediction horizon K on the average social distancing. The dotted cyan line refers to the number of currently available ICUs in Germany. The vertical dotted black line marks the end of social distancing measures for that setting.
leads to a reduction in the social distancing measures, as can be seen in Figure 14 .
In summary, increasing test capacities and/or ICU capacities helps to reduce measures like social distancing. However, the impact of the number of available ICUs appears to be much stronger. Nonetheless the qualitative shape of the solutions over time is not affected by varying these constraints.
In this paper, we demonstrated how mitigation of the COVID-19 epidemic can be achieved by a combination of age-stratified testing and social distancing measures while avoiding a breakdown of the health care system. We showed that in our compartmental model mass testing alone is insufficient to achieve this goal, as it would require unrealistic testing capacities.
As a remedy, we designed optimal social distancing strategies with a focus on applicability and acceptance in the general population, i.e. strategies with slowly changing contact reductions. The resulting social distancing measures imitate the measures actually taken in Germany, but are lifted at a much slower pace. Agedifferentiated contact reductions may improve upon these results as they yield qualitatively similar social distancing strategies and prioritize relaxing restrictions for the work-force and children.
To model the process of policy making more realistically, we used MPC which allows to adapt to deviations from the envisioned course of the epidemic by solving the optimal control problem repeatedly. Our analysis reveals that longer prediction horizons allow for faster lifting of restrictions although long-term predictions may be infeasible in practice. Additionally we showed that the amount of available intensive care units is a key factor influencing the required amount of social distancing.
We believe that our model with the chosen parameters reflects reality sufficiently well to provide qualitatively valid insight on how testing and social distancing can control the spread of SARS-CoV-2. We learned that mass testing alone is, assuming realistic testing capacities, not sufficient to avoid a breakdown of the health care system in Germany. To prevent this, one has to implement strict contact reductions early on, which, ideally, should then be eased slowly. If one allows these reductions to vary by age, one is able to relax restrictions for the (working) middle age group, at the cost of reducing contacts of the more vulnerable older population.
While short-term planning of measures is unable to control the exponential growth of cases, medium-term planning produces strategies that, qualitatively, do not differ from optimal ones while being flexible enough to adapt to new circumstances. Finally, as expected, the number of available intensive care units dictates how fast herd-immunity can be reached and how much total social distancing is necessary.
However, we caution the reader against interpreting these results in a quantitative way, as our model has not been devised to produce precise predictions. Similarly, we want to stress that we do not provide concrete policies to implement, as the impact of particular countermeasures on β is not easily quantified.
Concerning other influences on the epidemic's evolution, note that we have not yet considered vaccinations nor re-infections, both of which could be included in our model without difficulties, if parameters are available to model them. As our model is based on ODEs, interaction effects such as contact tracing cannot be included. Agent-based (stochastic) models are able to handle these critical effects and could be seen as a natural extension of our (deterministic) compartmental model. To solve the resulting stochastic optimal control problems would then require more sophisticated techniques, however.