Piecewise-constant optimal control strategies for controlling the outbreak of COVID-19 in the Irish population

We introduce a deterministic SEIR model and fit it to epidemiological data for the COVID-19 outbreak in Ireland. We couple the model to economic considerations — we formulate an optimal control problem in which the cost to the economy of the various non-pharmaceutical interventions is minimized, subject to hospital admissions never exceeding a threshold value corresponding to health-service capacity. Within the framework of the model, the optimal strategy of disease control is revealed to be one of disease suppression, rather than disease mitigation.


Introduction
Since late 2019, the COVID-19 respiratory disease originating from the SARS-CoV-2 virus [1] has caused a global pandemic, leading to over 9 million confirmed cases and over 660,000 confirmed deaths as of June 2020 [2]. As the epicentre of the pandemic has shifted from Asia to Europe and the Americas, governments across the world have taken steps to curb the spread of the disease. Due to the current lack of a vaccine for COVID-19, and the current very limited number of therapeutic agents available for treatment of COVID-19 [3], these efforts have taken the form of non-pharmaceutical interventions (NPIs).
Non-pharmaceutical interventions (NPIs) to curb the spread a pandemic like COVID-19 may take various forms, ranging from mass publicity campaigns to educate the public about proper hygiene, social distancing, testing, up to quarantining. The main purpose of such interventions is to reduce the effective reproduction number of the disease. Reducing the effective reproduction number serves either to mitigate the disease, or to suppress it entirely [4]. As per Ref. [4], we take 'suppression' to mean 'driving infections to the (unstable) = 0 fixed point', that is, driving the number of infectious cases to zero. Therefore, for the present purposes, 'suppression' and 'elimination' are synonymous. In contrast, mitigation serves to make the spread of the disease in the population more manageable (e.g. preventing a breakdown in health services due to an excess of hospital admissions). Mitigation may also involve a slow build-up of the disease in the population, with a view to achieving 'herd immunity' [5][6][7]. In contrast, suppression aims to reduce the effective reproduction number of the disease so far below unity that the disease dies out in the population.
In this work, we address the tradeoffs between mitigation and suppression for a specific country case (Ireland). We use official epidemiological data to generate a standard SEIR-type model for the * Corresponding author.
outbreak of COVID-19 in Ireland. Thereafter, we carry out a counterfactual exercise and look at what would be the optimal strategy of NPIs to reduce the effective reproduction number of the disease while simultaneously minimizing the cost to the national economy. The basic methodology is a standard compartmental SEIR model with multiple additional compartments for fitting the model to the extant data; the model is deterministic and based on coupled ordinary-differential equations. Thereafter, the model is reformulated in the context of Optimal Control Theory -the model is solved in such a way as to minimize a cost function, while satisfying certain state constraints which correspond to the requirement that hospital admissions remain below capacity for the duration of the outbreak.
Optimal control theory has already found great use in mathematical epidemiology. It has been used to compute optimal vaccination strategies for various diseases lethal to humans [8][9][10], as well as optimal containment strategies for animal diseases [11,12]. These works use the Pontryagin Maximum Principle [13] to formulate semi-analytical conditions for the optimal control strategy. In cases with a linear cost function (including optimal stopping time problems [11,12]), the optimal control is revealed to be a so-called bang-bang control, where the control switches between a maximum and a minimum value during at a finite number of switching points.
Our work is a departure from these prior works, in that we recognize the importance of state constraints in a context where it may be desirable to strictly control the number of infectious cases, with a view to preserving the integrity of a health system. Including state constraints is challenging from the mathematical point of view, and often a numerical approach is pursued [14]. Our work follows in this vein. In this work, we perform optimal control theory over a fixed time horizon of one year. The motivation is that this is a realistic time horizon over which to manage the spread of COVID-19 in population, before the advent of anticipated therapeutic interventions such as antiviral drugs and vaccines. This approach follows the recent work in Ref. [5]. This prior work is an important first step in the application of Optimal Control Theory to the COVID-19 outbreak. We build on this prior work by working with a linear cost function (Reference [5] uses a cost function that is quadratic in the intensity of the nonpharmaceutical interventions). A linear cost function may be more realistic from the point-of-view of economic modelling [15]. Ref. [5] involves modelling in the absence of state constraints; in the present article we include state constraints. This makes the mathematical analysis more involved, so a numerical optimization approach is favoured herein. Finally, we emphasize that although the model is calibrated for the Irish context, and the optimization performed in that context also, the methods are entirely generic, and will find use more broadly.
This article is organized as follows. In Section 2 we introduce the basic SEIR model and demonstrate how it is fit to epidemiological data for the COVID-19 outbreak in Ireland. In Section 3 we couple the model to economic considerations -we formulate an optimal control problem in which the cost to the economy of the various NPIs is minimized, subject to the satisfaction of certain key state constraints. Results are presented in Section 4, and conclusions are presented in Section 5.

Mathematical model
In this work, we develop an optimal control theory based on a deterministic compartmental model of epidemiology. Therefore, in this section we introduce an appropriate multi-compartment susceptibleexposed-infectious-removed (SEIR) model. The model is tailored to COVID-19, and is based on prior works [16,17]. Thus, the model includes specific compartments for susceptible individuals , exposed (but not yet infectious) individuals , and three separate infectious classes: infectious-pre-symptomatic ( ), infectious-asymptomatic ( ), and infectious-symptomatic ( ) individuals. Removed cases ( ) are considered also in a standard fashion. Finally, in order to fit the model to existing data, further compartments are introduced to the model to enable the recording of confirmed cases and deaths. In this way, we arrive at the following model. Time is measured in days.
Additional, non-standard compartments are introduced to allow for comparison with recorded data: Here, cases in the 'Ailing' compartment die and enter the -compartment at a rate . Also, it is assumed that all individuals entering the symptomatic compartment are tested -these cases enter the 'Awaiting Test' compartment, wait on average days for a test; results are subsequently reported and attributed to the -compartment. Including deaths in the model is important, as it enables us to fit the model to a wider range of epidemiological data involving both confirmed cases and deaths, thereby improving the fit of the model. Also, it opens up the possibility of performing optimal control on the model in such a way as to control or even limit the number of deaths over the course of the epidemic.
We further elaborate on the enumerated constants. The quantity is the total initial susceptible population. The coefficient denotes the average number of contacts between a susceptible person and an infectious person, per day. Correspondingly, the coefficient denotes the probability that an infected, pre-symptomatic person transmits the disease to a susceptible person during one of those contacts, and similarly for and . By allowing for different values between , , and , we can account for the possibility that asymptomatic cases are less infectious than symptomatic cases. We can also allow for the effect of case isolation, which will reduce the value of . The time (in days) denotes an incubation time when an exposed individual is not yet infectious, the time (in days) denotes the time interval during which an exposed individual is infectious before symptom onset. The coefficient denotes the fraction of pre-symptomatic cases that never show any symptoms of the disease. All infectious cases (whether asymptomatic or symptomatic) are assumed to remain infectious for a period days. It is assumed that all asymptomatic cases recover and are immune, while only 1 − of the symptomatic cases recover; the remaining become ailing and subsequently, die. Thus, represents the long-term case fatality rate of the disease.

Model fitting
For definiteness, we fit the model (1) to epidemiological data for the COVID-19 outbreak in Ireland. The epidemiological data is obtained from official sources [18], and involves both confirmed cases of COVID-19, and deaths attributed to COVID-19. To fit the data to the model (1), we have minimized the 2 norm, Here, labels time in days. Here also, = 76, corresponding to 77 days' worth of data, starting with the first recorded case of COVID-19 in Ireland on February 29th 2020 ( = 0) and ending on May 17th 2020 ( = 76). The justification for choosing this time-period is that it includes the initial phase of the COVID-19 outbreak before diseasecontrol measures were introduced (beginning 12th March 2020). The chosen time-period ends before the beginning of the phased lifting of these controls (May 18th 2020). Finally, 1 denotes the day with the first recorded death and 2 denotes the first day in the outbreak after which the number of new recorded deaths per day is strictly greater than zero.
We use simulated annealing to minimize as a function of the model parameters. To reduce the complexity of the problem slightly, we re-write Eq. (1a) as: a similar modification applies to Eq. (1b). Here, 0 = , = ∕ , and = ∕ . As such, and should henceforth be understood as Table 1 Effect of various interventions on the value of the basic reproductive number 0 . The formula for 0 is given in Eq. (8), with 0 replaced by . The given values of correspond to 1 and 2 respectively in the optimization problem (6), the numerical values of which are given in Appendix A. Similarly, the value for the case-isolation parameter is obtained from the optimization problem (6), with the details provided again in Appendix A. ratios of probabilities -e.g. should be taken to mean the probability that an asymptomatic individual transmits the disease to a susceptible individual, divided by the probability that a pre-symptomatic individual transmits the disease to a susceptible individual. A similar interpretation should henceforth be applied to .
At the same time, the optimization problem is rendered more complex, in several respects. First, the model is initialized at = − , with and all other compartments equal to zero. Here, is an unknown parameter, In this context, the total population is assumed to be susceptible to the disease, and ( = − ) = 1 corresponds to 'patient zero'. Furthermore, 0 changes over the course of the epidemic, corresponding to the implementation of disease controls [19]. To account for this, we replace 0 in the model (1) by: In the first time period ( < 13), no disease controls are implemented. In the second time period (13 ≤ < 28), all bars, schools and universities are closed, and mass gatherings are banned. Finally, in the third time period ≥ 28, all non-essential services and industries are closed down, and travel is restricted to within 2 km of a person's house, with certain limited exceptions. Eq. (5) already contains some simplifications that do not correspond exactly to the actual course of the COVID-19 control efforts in Ireland -e.g. a four-day period (12th March 2020-15th March 2020, inclusive) where schools and universities were shut down and mass gatherings banned but bars remained open. Also, after May 5th 2020 ( = 66), the 2 km restriction was extended to 5 km. However, we ignore these details, as the main purpose of this article is not necessarily to track and predict the course of COVID-19 in Ireland, but rather to develop a simple, generic methodology to compute the optimum application of non-pharmaceutical interventions to control the spread of COVID-19. As such, we seek to solve the optimization problem for which is a 13-dimensional optimization problem. We solve the optimization problem (6) using simulated annealing. Bounds for the search space are estimated from the literature -these are reported in Appendix A, where we also provide the optimal parameter values corresponding to Eq. (6). Results are shown in Fig. 1. Qualitative agreement between the SEIR model (1) and the data can be seen in the figure. To demonstrate quantitatively the agreement between the model, the data, and prior characterizations of the epidemiology of COVID-19, we have computed the basic reproduction number of the model (1) using the maximum eigenvalue of the nextgeneration matrix [20], 0 = max spec( −1 ). For the model (1), the matrices and are given by: hence (the reason for including the explicit functional dependence on the parameters 0 and will become clear in what follows). We have used the parameters obtained from solving the optimization problem (6) to compute the value of 0 . With case isolation (i.e. = 0.3405, as per the solution of Eq. (6)), we have computed 0 = 2.6945. Without case isolation (i.e. = 1), we have computed 0 = 3.4495. These results are consistent with an existing estimates of the basic reproduction number for COVID-19 [21]. Furthermore, by replacing 0 in Eq. (8) with both 1 and 2 , we can estimate the effect of various measures to control the spread of the epidemic. This is done in Table 1 -the results are very close to other independent estimates of the same quantities [19]. These findings confirm the robustness of the basic mathematical model (1) and therefore provide a solid foundation on which to basic a theory of optimal control in the remaining sections.
Finally, we emphasize that the focus of this article is on developing a numerical framework for optimal control for epidemic outbreaks. Therefore, quantifying the uncertainty in the fit of the SEIR model (1) to the epidemiological data takes a back seat. However, in Appendix A, along with the optimal parameter values corresponding to Eq. (6), we also provide the confidence intervals for the same, as well as confidence intervals for the basic reproductive number 0 and the effective reproductive numbers in Table 1.

Optimal control theory
In this Section we model the imposition of non-pharmaceutical interventions (NPIs) by replacing the constant value 0 in Eq. (1) with a piecewise-constant function denoted by ( ). We solve the equations from the initial time = − (initial conditions (4)) out to a fixed final time = . Specifically, we take: Here, ( ) is a second piecewise-constant functions. In this way, the model (1) is solved with no controls for ≤ 0 ; controls are introduced at time 0 . The time 0 is fixed. For consistency, the following bounds apply to ( ): The piecewise-constant function ( ) is characterized as follows: Here, is an integer (one or greater), 1 , … , are real numbers between zero and one, and 1 , … , , −1 are switching times (real numbers between 0 , such that 1 < ⋯ , −1 < ). In practice, the switching times correspond to how often the NPIs are updated. We recall that 0 = , where is the average number of contacts between a susceptible person and an infections person, per day, and denotes the probability that an infected, pre-symptomatic person transmits the disease to a susceptible person during one of those contacts. Therefore, the function ( ) describes at a population level how to implement a strategy of NPIs -through a reduction of either or . For example, ( ) = 0.5 may correspond to a 50% reduction in , resulting from a corresponding reduction in mobility for the susceptible population.

Formulation
We quantify the cost to the economy of implementing a sequence { 1 , … , } of controls by the cost function with the convention that ,0 = 0 and , = . As such, is a function of the real variables { 1 , … , , 1 , … , , −1 }. The rationale for this choice is as follows. The aim of each NPI is to reduce the transmission of the virus. Other than the 'easy' interventions such as public information campaigns, handwashing, etc. the interventions aim to reduce the number of daily contacts between individuals. Assuming economic activity is proportional to this number of contacts, there is a natural linear relationship between the level of the NPIs and the cost to the economy (e.g. if people's daily contracts are reduced by 10%, then one can expect economic activity to be reduced by the same factor). However, we recognize that some of the more intense NPIs may carry a disproportionate cost to the economy and for that reason, the linear relationship in Eq. (11) may break down as → 1. For that reason, in what follows, we investigate the robustness of our results to the assumption of a linear cost function (Section 4).
We furthermore insist on an optimal control problem where human life is put on a very high footing. As such, we propose the following additional constraints on the optimal control problem: In this context, the positive constant may be thought of as the percentage of symptomatic cases who require a hospital bed (or a bed in ICU) at time , and the positive constant represents a corresponding capacity limit (number of hospital beds, number of ICU beds, etc.). Finally, we impose the additional constraint which is true if and only if −1 (1 − ) − −1 ≤ 0 at = . This rules out any 'optimal' strategy in which a large epidemic peak would occur just beyond the horizon at = . The exclusion of such strategies is desirable for public-health reasons. As such, the optimal control problem to be solved is the following: Under the epidemic modelled by Eq. (1), over a time horizon , compute the minimizer of subject to the constraints (12)-(13).

Discussion
We discuss the optimal control problem (14) in a more general context. We use the notation = ( , , , , , ). We look at the optimal control problem Subject to d ∕d satisfies Eq. (1), 0 ≤ ≤ , The control is activated at = 0; between = − and = 0, the evolution of d ∕d is given by Eq. (1) with = 0. This particular optimal control problem can be put into the framework of Theorem 23.11 in Ref. [22] -this framework provides the conditions such that, if there is at least one admissible process ( , ) for which is finite, the optimal control problem admits a solution. The structure of the problem (15) is relatively simple, and to show that the conditions of Theorem 23.11 in Ref. [22] apply here, it suffices to show that solutions ( ) of the ordinary differential equation (1) remain in a bounded hypercube in R 6 ; this is done in Appendix C. Furthermore, if the number of infectious individuals is sufficiently small at = 0, then by making ( ) = Const. sufficiently large, the solution ( ) will follow a standard epidemic curve with ( ) ≤ -this gives the required admissible process. This framework gives very general necessary conditions for a solution to the problem (15) to exist. Often, standard analytical approaches can be applied to compute this solution analytically. The application of these approaches is hampered in the present context by the presence of the state constraints in Eqs. (12)-(13). For a system of ordinary differential equations ∕ = ( , ), state constraints are generically expressed as ( ) ≤ 0 (for this brief presentation we assume for discussion purposes that there is only one state constraint, hence is a scalar-valued function of its arguments). The order of the constraint is the minimum value of for which ∕ can be expressed as a linear function of , that is, Here, ∕ = (∇ ) ⋅ , and and are continuous functions of their arguments. In the absence of such state constraints, a convenient way of solving optimal control problems is the Pontryagin Maximum Principle [13]. For cost functions that are linear in the control ( ), a consequence of the Pontryagin Maximum Principle is that the optimal control is bangbang. That is, if 0 ≤ ( ) ≤ , then the optimal control switches between ( ) = 0 and ( ) = , with a finite number of switching times. With state constraints, the optimal control tends to be a mixture of bang-bang control sequences and 'boundary arcs', during which time the control ( ) is such that the system maintains ∕ = 0 [23]. The state constraints in Problem (14) are of high order (for instance, Eq. (12) is a third-order state constraint). The analysis of such problems is demanding, as the constraints involve the analysis of new multipliers associated with each individual constraint [14]. For this reason, the focus of this work is on computing optimal controls numerically. Clearly, the resulting absence of theoretical analysis means that our  results should be treated with care since there is no guarantee that they are optimal. However, we have carefully examined the robustness of our numerical calculations -this is described in detail in Section 4 and also, in Appendix B.
There is also a second and more compelling reason for our use of numerical optimization methods in this work. From the outset we rule out the search for optimal controls that involve boundary arcs. Although this restriction may result in the ruling out of globally optimum solutions, we are guided by the intended application of the work, where the outbreak of an epidemic is to be controlled using public-health interventions. These interventions tend to be on the population level. Communication with, and collaboration from, the entire population is necessary for such measures to work. As such, piecewise-controls of the type (10) are more realistic. The formulation of the optimization problem in terms of such discrete controls lends itself naturally to a numerical approach.

Methodology
We solve the optimization problem (14) using simulated annealing. We use the built-in simulated annealing algorithms in Matlab/Octave. The source code is obtainable via an online repository [24]. During each round of simulated annealing, a random set of controls { 1 , … , , 1 , … , , −1 } is drawn from the feasible values. The model (1) is then solved numerically using ODE45 in Matlab/Octave and the penalty function (11) is calculated. Extra terms are added to the penalty function for cases where the choice of controls causes the constraints (12)-(13) to be breached. The simulated annealing algorithm generates successive random sets of controls until a global minimum value of the penalty function is attained. Simulated annealing is guaranteed asymptotically to converge to the global minimummore precisely, for any finite-dimensional problem, the probability that the simulated annealing algorithm terminates with a global optimal solution approaches 1 as the annealing schedule is extended [25]. In practice, the convergence can be very slow, and the Matlab/Octave codes employ built-in stopping criteria to confirm if a global minimum has been reached. These built-in stopping criteria can occasionally cause the simulated-annealing algorithm to converge before the attainment of the global minimum. We have checked this by comparing the results of the simulated annealing with other fundamentally different optimization methods with independent stopping criteria (e.g. particleswarm optimization), and the results are the same. This inspires confidence in our numerical method, although the absence of theoretical analysis of the optimal control problem (14) at this time means that our results should be treated with some care.

Numerical results
We solve a variant of optimization problem (14) numerically. Specifically, for various values of the number of controls , we compute subject to the constraints (12)-(13) (we also compute the corresponding values { 1 , … , } and { 1 , … , , −1 } which achieve the minimum).
We report on the trends as the integer is increased from = 2 to = 12. The model (1) is solved subject to the initial conditions (4), with = 0 corresponding to February 29th 2020, the day on which the first case of COVID-19 was recorded in Ireland. The model is solved without any controls until = 0 = 13 Days (cf. Eq. (5)), whereupon the control ( ) is activated. In this way we look at counter-factual scenarios, to determine what would have been the optimal sequence of controls to implement, to control the outbreak of COVID-19 in Ireland. We fix the final time to be = 0 + 365 days. In this, we are motivated by similar recent work which looks at how to control an outbreak of COVID-19 over a timeframe prior to the development of pharmaceutical treatments of the disease [5]. Finally, based on the official statistics [26], we use = 0.016 and = 300 in the state constraints (12)- (13), this corresponds to an ICU admission rate of 1.6% among all confirmed cases of COVID-19.
To take account of strategies which reduce the spread of the disease to a very low level, we supplement Eq. (3) by introducing an indicator function, Then, Eq. (3) is further modified to read, = −( ) 0 ( + + ); a similar modification applies to the equation for the 'exposed' compartment. In this way, once the disease is reduced to a low level in the population, in the sense that + + < 1, an unrealistic resurgence in the disease by a notional 'fraction of a person' is ruled out.
The main purpose of this section is to characterize the solution of Eq. (16) as a function of , where  That is, we look at the structure of the solution of the problem (16) for various values of . We deliberately restrict ourselves to ≤ 0.8. When the model (1) was fitted to the data for Ireland, ≈ 0.8 was the maximum reduction in achieved, under a severe set of interventions designed to curb the spread of COVID-19 (all non-essential services and industries closed down, travel restricted to within 2 km of a person's house, with certain exceptions).

One switching time, = 2
We start by looking at the structure of the solution of the problem (16), for a range of values of , and for = 2. Although switching between different sets of controls only once in the time period might be unrealistic, this is a first step that helps with understanding (we investigate > 2 in what follows). As such, we plot the minimum value as a function of in Fig. 2. For large values of , 0.7, the optimal control is bang-bang with one switch: the control is enabled at the maximum value ( ) = for between 58 ( = 0.8) and 193 days ( = 0.71). The evolution of the infection under these controls is shown in Fig. 3 -under these bangbang controls, the epidemic is seen to be suppressed, and the number of infectious cases is zero at the end-time . In contrast, for 0.7, the optimal control is no longer bang-bang. Instead, the optimal control strategy is revealed to be independent of : for an initial time  interval, the value ≈ 0.607 is optimal, during which the boundary of the state constraint ( ) = is attained, for some . In a second time interval, the value ≈ 0.468 is optimal, during which time the boundary of the state constraint is again attained. The switching time occurs at 1 − 0 ≈ 256.5 Days. The evolution of the epidemic under this scenario is shown in Fig. 4. This scenario corresponds to a 'containment' of the epidemic: the number of infectious cases is non-zero at the end-time , but the maximum capacity of the hospital system is not breached. Finally, a comparison between a disease outbreak with controls, and one without any controls, is shown in Fig. 5. The imposition of controls contains or suppresses the epidemic: specifically, max( ) ≈ 3.7 × 10 5 with no controls, max( ) = ∕ = 18,750 with = 0.7, and max( ) ≈ 100 with the optimal control = 0.8. Comparison between these scenarios and the actual course of the outbreak is shown in Fig. 6. The model predictions for the counterfactual scenario of an unmitigated outbreak are consistent with the initial modelling of the epidemiology of COVID-19 elsewhere, e.g. Ref. [4].
In the case with one switching time (i.e. = 2), the optimization problem (16) contains only three parameters, and the optimal control strategies found in Figs. 3-4 can be investigated separately via a 'brute-force approach'. This is desirable because it helps with our understanding and also, because it gives further confidence in the correctness of the simulated-annealing approach in isolating the global minimum. As such, we have evaluated for a wide range of values, by direct computation. By fixing 1 , a two-dimensional plot can be made, and a value ( 1 ) = min ] (subject to the state constraints) can be obtained graphically. By repeating this procedure over all values of 1 , the optimal control strategies reported in Figs. 3-4 are recovered. Sample two-dimensional plots along these lines are shown in Fig. 7. The minimum marked in panel (a) corresponds to the global minimum with the imposed constraint = 0.8 on the control. This corresponds to a bang-bang control with = with one switching time at 1 − 0 ≈ 0.198 × 365 Days, and = 0 thereafter. It can be asked if the minimum ( 1 = 0.8, 2 = 0, 1 − 0 ≈ 0.198 × 365) is a true bangbang control or is instead a boundary control, as this optimal control appears to occur along the curve max ∈[ 0 , ] ( ) = in Fig. 7(a). The answer is that the surface ( 1 , 2 ) = max ∈[ 0 , ] ( ) − possesses a jump discontinuity, which arises as a result of the extra condition (17)- (18) in the model (1). Therefore, the optimal control ( 1 = 0.8, 2 = 0, 1 − 0 ≈ 0.198 × 365 Days) occurs on one side of a jump discontinuity and thus produces a solution that never attains the boundary value max ∈[ 0 , ] ( ) = . This explanation can be verified by removing the condition (17)-(18) from the model; in that case, the discontinuity goes away, and the optimal control ( 1 = 0.8, 2 = 0, 1 − 0 ≈ 0.198×365 Days) shifts elsewhere. However, we rule out further study of this special case as the condition (17)-(18) is deemed necessary for model realism.
We emphasize finally that although the two control strategies represented by Figs. 3-4 each represent a solution to the optimal control  problem (16), they are by no means equivalent. The 'bang-bang' control strategy (high ) comes with a much lower economic cost (lower ). Therefore, subject to the model assumptions and the framework of the optimal control problem (14), an eradication strategy is far preferable to a containment strategy.

Several switching times, > 2
We extend the foregoing analysis to include control strategies with multiple switching times. The results are summarized succinctly in Fig. 8. The results are qualitatively similar to what is observed in the case of a single-switch ( = 2) control -that is, for large values of , the bang-bang control with a single switch and the complete suppression of the disease is optimal, while for lower values of a containment strategy is optimal. Under the mitigation strategy, the state constraint max ∈[ 0 , ] ( ) = is attained at several times -the attainment of the bound corresponds to local maxima in the function ( ) of infectious cases (e.g. Fig. 9). Also, the cost to the economy of the mitigation strategy is reduced as the number of switching times is increased. Overall, however, the strategy of complete suppression is optimal: the higher the value of permitted, the lower the economic cost of the control measures.
The results in Fig. 8 are for ≤ 12; we do not rule out that for > 12, a mitigation strategy may become competitive with a suppression strategy. However, all mitigation strategies (no matter how many switches) have an inherent instability-like behaviour -small deviations from the optimal control lead to large deviations in the outcome of the epidemic. This is illustrated in Fig. 10, where we show 10 different simulation results corresponding to control sequences that are within 10% of the optimal control. In several cases, this small deviation away from the optimal control leads to a breach of the state constraint concerning ICU capacity by a factor of 4. Thus, the maintenance of the state constraints is highly sensitive to small deviations away from the optimal control. Now, the non-pharmaceutical interventions discussed in this article are population-based interventions, concerning such public-health initiatives as social distancing, testing suspect cases, quarantining infected individuals, up to so-called 'lockdown' of entire populations. It can be appreciated that these are blunt tools for epidemic control, and are dependent on a high degree of compliance from the public. Consequently, a deviation of 10% in the implementation of these controls seems reasonable. And yet such a small deviation leads to a four-fold breach in the ICU capacity. This sensitivity in the outcomes to small deviations in the controls is another reason why the mitigation strategy is inferior to the suppression strategy.

Quadratic penalty function
We have repeated the analysis of the previous sections for a quadratic cost function, The motivation for considering Eq. (21) is to look at the hypothetical scenario where the cost to the economy of the non-pharmaceutical interventions is nonlinear in the intensity of the interventions. In this way, Eq. (21) describes a scenario where the more intense interventions ( = 1) carry a disproportionately severe cost in comparison to less intense interventions ( < 1). In spite of the differences between the linear and quadratic cost functions, the results in both cases are qualitatively the same -e.g. a plot of the minimum value of as a function of for the quadratic case is shown in Fig. 11, to be compared with the linear cost function in Fig. 8. The main difference between these two figures is that the crossover where the optimal strategy switches from mitigation to suppression occurs at a higher value of for the quadratic cost function -that is, under the assumption of the quadratic cost function, mitigation is an optimal strategy for a wider range of values of . That being said, the sensitivity of the mitigation strategy to small deviations in the controls still applies for the case of the quadratic penalty function.

Sensitivity of results to parameter variations
We briefly discuss the sensitivity of our results to parameter variations, focusing on some key parameters of interest. In the first instance, we remark that we have taken 0 as fixed in this work -that is, the day on which NPIs are first introduced is fixed a priori. We have taken ( ) under an optimal-control strategy in which the boundary is attained. Control parameters: = 6, = 0.68. Other curves: plot of ( ) for control strategies that are within 10% of the optimal control strategyspecifically, each and is varied independently by 10% with respect to the optimal values. Fig. 11. Plot of as a function of , with ≥ 2, quadratic cost function (Eq. (21)).
In the same way, we emphasize again (cf. Section 1) that the time horizon of the optimal control problem is fixed, with = 1 year. The motivation is that this is a realistic time horizon over which to manage the spread of COVID-19 in population, before the advent of anticipated therapeutic interventions such as anti-viral drugs and vaccines. Introducing a variable time horizon for the optimal control problem can be addressed in future work -the same observations as above will apply, namely the importance of suppressing the initial exponential growth of the epidemic.
Finally, it can be noted that optimal control solutions are well known to suffer from low robustness. The same caveat applies here -it may be anticipated that the control parameters { 1 , … , , 1 , … , , −1 } will vary under small changes in the model parameters (for instance, parameter variations within the bounds of the confidence intervals in Appendix A). As such, the optimal control parameters are not expected to be prescriptive. Instead, the general trends in the results are more important -in particular, the importance of suppressing any initial sharp rise in cases, such that high-intensity NPIs can be lifted as soon as possible, with a view to minimizing the economic cost of the suppression strategy.

Conclusions
Summarizing, we have introduced a mathematical model to describe the outbreak of COVID-19 in Ireland in March-May 2020. Based on this model, which we compare to actual epidemiological data, we have used numerical methods from optimal control theory to answer a hypothetical question -what would have been the optimal strategy for the control of the epidemic? The optimal control theory is based on minimizing a cost function which is related to the cost to the economy of the various control measures. At the same time, the theory introduces state constraints as a means of maintaining the number of infectious cases below a maximum threshold such that hospital capacity is not overwhelmed. By insisting on piecewise-constant controls (appropriate for modelling rather crude population-level public-health interventions) and state constraints, the optimal-control model is not amenable to theoretical analysis -hence a numerical approach is considered. As such, we have performed the necessary optimizations via numerical simulated annealing -the results show that a suppression strategy (implemented promptly) is superior to a mitigation strategy.
The optimal suppression strategy may be compared to the strategy actually implemented by the Irish Government in the early stages of the epidemic (March-May 2020). Broadly, the optimal strategy is similar to the actually implemented strategy, in the sense that both involve high-intensity NPIs at the beginning of the epidemic, with a view to suppressing the outbreak. However, the optimal strategy envisages maximum-intensity NPIs at the very outset -this saves lives (e.g. Fig. 6, with 1547 deaths at − 0 = 65 of the outbreak in real life, and 57 deaths on the same day under the optimal strategy, and an asymptotic final number of deaths of 61 under the same optimal strategy). Therefore, imposing high-intensity NPIs from the very beginning saves lives; it is also economically less costly, under the penalty function formulated in this work.
The optimal strategy predicted by the model envisages the maximum level of non-pharmaceutical interventions (NPIs) from the outset (Day 13 in our calculations, cf. Eq. (5)), and lasting 58 days at = 0.8. In contrast, the strategy implemented by the Irish Government involved some controls on Day 13, maximum controls on Day 28, followed by gradual easing of controls on Day 79 (May 18th 2020). The period between the initial NPIs and the first easing of NPIs is therefore 79−13 = 66 Days. Thus, the immediate imposition of maximum controls on Day 13 would have resulted in the lifting of the most intense NPIs roughly 14 days sooner. Naturally, this estimate is subject to model uncertainty, but the clear message (supported by the other findings in this paper) is that prompt implementation of maximum NPIs is best.
Finally, it can be emphasized that the coupling of the epidemiological model to the economy in the present model framework is rather crude (e.g. the linear cost function (11)). In future, it can be envisaged that epidemiological models can may be coupled to economic models of production [15] or employment [27] to better understand the interplay between the different societal costs of a disease outbreak. Equally, the approach presented in this work may form a basis for the design of an optimal vaccination strategy if and when a vaccine becomes available for COVID-19.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Optimal value is checked to make sure it is strictly within the bounds 0 1 ≈ due to occurrence of peak viral load close to symptom onset [31]. Also, < [28]. Hence, = ∕ ≤ 1. In this Appendix we provide the necessary background information to fit the model (1) to the epidemiological data in Ref. [18] for the outbreak of COVID-19 in Ireland. A nonlinear optimization is performed to minimize the penalty function ( , 0 , 1 , 2 , , , , , , , , ) in Eq. (6). First, bounds on the range of acceptable parameter values are obtained in Table A.2. The resulting optimal parameter values are provided in Table A.3. In this table, we also present confidence intervals for these optimal parameter values -these have been obtained by bootstrapping -the model residuals are resampled with replacement times. These residuals are then used to compute synthetic data sets, and the model (1) is fitted to each data set in turn. This then gives a range of 'optimal' parameter values; a histogram based on this range of values is constructed and the 2.5 and 97.5 percentiles are computed to produce the reported confidence intervals. We use = 100 and = 1000, with little or no difference observed between these two choices.
Furthermore, we have used the data in Table A.3 to generate estimates of the basic reproduction number and the effective reproduction number under various NPIs, with appropriate lower and upper bounds. The results are shown in Table A.4, these may be compared with Table 1 in the main paper.
According to the estimates in Table A.2, the mean value of the infectious periods + is between 2.15 and 5.45 days. This is consistent with Ref. [16]. However, some other studies indicate that the infectious period may be significantly longer (e.g. up to 10 days [31]). As the main purpose of this article is to establish a methodology for computing optimal control strategies, this discrepancy is not important in the current context. However, exploring the effect of this range of values on the results may be of interest in future work.

Appendix B. Description of simulated annealing method for optimization
To maximize the chance of finding a globally optimal control, we use the method of Simulated Annealing. We exploit the builtin simulated annealing capabilities in Matlab and Octave for these purposes. As such, we solve Eq. (1) for an arbitrary sequence of controls { 1 , … , , 1 , … , , −1 }. In Matlab/Octave programming we execute the following command:

penalty=ode_solve_seir(u);
Here, ode_solve_seir is a Matlab/Octave ODE45 solver which takes in the input = { 1 , … , , 1 , … , , −1 }, solves Eqs. (1) out the final time = , and returns the penalty function (11). The state constraints (12)- (13)  We use the default convergence criteria in the simulated annealing algorithm in Matlab. In particular, the algorithm is deemed to have converged if the average change in the objective function over the previous 500 × (Problem Size) is less than 10 −6 . Here, (Problem Size) is the size of the optimization problem which in the present case is 2 − 1.
We have systematically varied these parameters and find no significant change to the estimated global minima. In order to carry out large-scale simulated annealing runs (e.g. over 40 runs to generate the graph in Fig. 8, some taking over 3 days), we have implemented the simulated annealing algorithms in Octave (a GNU open-source alternative to Matlab) and carried out largescale runs on a cluster. The cluster is made up of 10-core/20-thread Intel(R) 'Ivy Bridge' Xeon(R) CPUs (2.50 GHz, 32 GB RAM). Each optimization job is run on a single core; however, the jobs are carried out simultaneously in batch mode. The default convergence criteria for the simulated annealing algorithm are different from those in Matlab [32] (the average change in the objective function over the 5 previous pseudo-temperature changes in the simulated annealing algorithm should be less than 10 −10 ). The results under the Matlab convergence criteria and the Octave convergence criteria are near-identical (exactly-identical results are unlikely, since simulated annealing is a non-deterministic algorithm). This inspires confidence in the robustness of the results generated in the paper, although certainly, without a theoretical analysis of the optimal control problem, convergence to a global minimum is never guaranteed.

Appendix C. Solutions of the SEIR model remain in a bounded set
We show that solutions ( , , , , , ) of Eq. (1) remain bounded for all time, once suitable initial conditions are provided. For clarity of the exposition, we assume that Eq. (1) is posed on a time interval (0, ], with initial conditions provided at = 0; specifically, Since ( + + + + + ) ( ) = , it follows that each of ( , , , , , ) remains between 0 and for all ∈ [0, ] and hence, the solution of the ordinary differential equation (1) is bounded.