key: cord-0197654-rkprvqey
authors: Liu, Hailiang; Tian, Xuping
title: Data-driven optimal control of a SEIR model for COVID-19
date: 2020-12-01
journal: nan
DOI: nan
sha: 3878019d52d2f4f163d1749ed41ec72177c8e86d
doc_id: 197654
cord_uid: rkprvqey

Since the first case of COVID-19 in December 2019, a total of $11,357,322$ confirmed cases in the US and $55,624,562$ confirmed cases worldwide have been reported, up to November 17, 2020, evoking fear locally and internationally. In particular, the coronavirus is surging nationwide at this time. In this paper, we investigate a basic Susceptible-Exposed-Infectious-Recovered (SEIR) model and calibrate the model parameters to the reported data. We also attempt to forecast the evolution of the outbreak over a relatively short time period and provide scheduled optimal controls of the COVID-19 epidemic. We provide efficient numerical algorithms based on a generalized Pontryagin Maximum Principle associated with the optimal control theory. Numerical experiments demonstrate the effective performance of the proposed model and its numerical approximations.

The outbreak of COVID-19 epidemic has resulted in over millions of confirmed and death cases, evoking fear locally and internationally. It has a huge impact on global economy as well as everyone's daily life. Numerous mathematical models are being produced to forecast the future of COVID-19 epidemics in the US and worldwide. These predictions have far-reaching consequences regarding how quickly and how strongly governments move to curb the epidemic. Mathematical models can be profoundly helpful tools to make public health decisions and ensure optimal use of resources to reduce the morbidity and mortality associated with the COVID-19 pandemic, but only if they are rigorously evaluated and validated and their projections are robust and reliable.

One of the well-known and basic models in epidemiology is the SIR model proposed by Kermack and McKendrick [17] in 1927. Here, S, I, R represent the number of susceptible, infected and recovered people respectively. They use an ODE system to describe the transmission dynamics of infectious diseases among the population. In the current COVID-19 pandemic, actions such as travel restrictions, physical distancing and self-quarantine are taken to slow down the spread of the epidemic. Typically, there is a significant incubation period during which individuals have been infected but are not yet infectious themselves. During this period the individual is in compartment E (for exposed), the resulting models are of SEIR or SEIRS type, respectively, depending on whether the acquired immunity is permanent or otherwise. Also such models can show how different public health interventions may affect the outcome of the epidemic, and can also predict future growth patterns.

In this paper, we integrate the optimal control with a specific SEIR model, though the developed methods can be readily adapted to other epidemic ODE models. More precisely, we introduce a dynamic control model for monitoring the virus propagation. Here the goal is to advance our understanding of virus propagation by learning the model parameters so that the error between the reported cases and the solution to the SEIR model is minimized. In short, we formulate the following optimization problem

L(U (t i )) + g(U (T )), s.t.U = F (U, θ) t ∈ (0, T ], U (0) = U 0 , θ ∈ Θ. HereU = d dt U (t), θ is a time-varying vector of model parameters, Θ is the admissible set for model parameters θ, and U = [S, E, I, R, D] corresponds to the susceptible, exposed, infected, recovered and deceased population. The loss function J is composed of two parts: L measures the error between the candidate solution to the SEIR system and the reported data at intermediate observation times, and g measures the error between the candidate solution and the scheduled control data at the end time. This is in contrast to the standard approach of fitting the model to the epidemic data through the non-linear least squares approach.

In the present work we provide a fairly complete modeling discussion on parameter learning, prediction and control of epidemics spread based on the SEIR model. We begin with our discussion of basic properties of the SEIR model. Then we show how the parameters can be updated recursively by gradient-based methods for a short time period while parameters are near constants, where the loss gradient can be obtained by using both state and co-state dynamics. For an extended time period with reported data available at intermediate observation times, the parameters are typically time-varying. In this general setting, we derive necessary conditions to achieve optimal control in terms of the chosen objectives. The conditions are essentially the classical Pontryagin maximum principal (PMP) [29] . The main differences are in the way we apply the principle in each time interval, and connect them consistently by re-setting the co-state V at the end of each time interval. We thus named it the generalized PMP. We further present an algorithm to find the numerical solution to the generalized PMP in the spirit of the method of successive approximations (MSA) [9] . Our algorithm is mainly fulfilled by three parts: (1) discretization of the forward problem in such a way that solutions to U remain positive for an arbitrarily step size h; (2) discretization of the co-state equation for V is made unconditionally stable; (3) Minimiziation of the Hamiltonian is given explicitly based on the structure of the SEIR model. This data-driven optimal control approach can be applied to other epidemic models. In particular, the prediction and control can be combined into one framework. To this end, the cost includes terms measuring the error between confirmed cases (infection and death) and those predicted from the model during the evolution, and terms measuring the error between the scheduled numbers as desired and those predicted from the model during the control period.

1.2. Further related work. In the mathematical study of SIR models there is an interplay between the dynamics of the disease and that of the total population, see [2, 27, 12, 33] . We refer the reader to [25, 11] for references on SEIRS models with constant total population and to [24] for the proof of the global stability of a unique endemic equilibrium of a SEIR model. Global stability of the endemic equilibrium for the SEIR model with non-constant population is more subtle, see [23] . Epidemic models have been treated using optimal control theory, with major control measures on medicine (vaccination), see e.g., [31, 19, 22, 15, 21, 28, 3] , and for Cancer immunotherapy control [7] .

For COVID-19, a data-driven model has been proposed and simulated in [26] , where both the SIR model and the feed-forward network are trained jointly. Compared to this work, our data-driven model is to consider the optimal control for SEIR with rigorous derivation of optimality conditions and stable numerical approximations. On the other hand, our data-driven optimal control algorithm may be interpreted as training a deep neural network in which the SEIR model serves as the neural network with parameters to be learned. This is in contrast to the study of residual neural networks using neural ODEs (see e.g., [8] ). SIR models with spatial effects have also been studied [4, 20] . A SIRT model was proposed in [4] to study the effects of the presence of a road on the spatial propagation of the COVID-19. Introduced in [20] is a mean-field game algorithm for a spatial SIR model in controlling propagation of epidemics.

1.3. Organization. In Section 2, we motivate and present the SEIR system for modelling the time evolution of epidemics, and integrate it with data into an optimal control system. Both algorithms and related numerical approximations are detailed in Section 3. Section 4 provides experimental results to show the performance of our data-driven optimal control algorithms. We end with concluding remarks in Section 5.

2. The SEIR model and optimal control 2.1. Model formulation. A population of size N (t) is partitioned into subclasses of individuals who are susceptible, exposed (infected but not yet infectious), infectious, and recovered, with sizes denoted by S(t), E(t), I(t), R(t), respectively. The sum E + I is the total infected population.

The dynamical transfer of the population is governed by an ODE systeṁ

subject to initial conditions S 0 , E 0 , I 0 , R 0 , D 0 . Here D(t) denotes deaths at time t, A denotes the recruitment rate of the population due to birth and immigration. It is assumed that all newborns are susceptible and vertical transmission can be neglected. d is the natural death rate, µ is the rate for virus-related death. γ is the rate for recovery with 1/γ being the mean infectious period, and is the rate at which the exposed individuals become infectious with 1/ being the mean incubation period. The recovered individuals are assumed to acquire permanent immunity (yet to be further confirmed for COVID-19); there is no transfer from the R class back to the S class. β is the effective contact rate. In the limit when → ∞, the SEIR model becomes a SIR model. Note that the involved parameters do not correspond to the actual per day recovery and mortality rates as the new cases of recovered and deaths come from infected cases several days back in time. However, one can attempt to provide some coarse estimations of the "effective/apparent" values of these epidemiological parameters based on the reported confirmed cases.

Solution properties of the SEIR model. Classical disease transmission models typically have at most one endemic equilibrium. If there is no endemic equilibrium, diseases will disappear. Otherwise, the disease will be persistent irrespective of initial positions. Large outbreaks tend to the persistence of an endemic state and small outbreaks tend to the extinction of the diseases.

To understand solution properties of the SEIR system, we simply take A = bN with b the natural birth rate, and focus on the following sub-systeṁ

with the total population size N = S + E + I + R. By adding the equations above we obtainṄ

Let s = S/N, e = E/N, i = I/N, r = R/N denote the fractions of the classes S, E, I, R in the population, respectively, then one can verify thaṫ

From biological considerations, we study system (2.2) in a feasible region

It can be verified that Σ is positively invariant with respect to the underlying dynamic system. We denote ∂Σ and Σ 0 the boundary and the interior of Σ in R 3 , respectively. A special solution of form P 0 = (1, 0, 0) on the boundary of Σ is the disease-free equilibrium of system (2.2) and it exists for all non-negative values of its parameters. Any equilibrium in the interior Σ 0 of Σ corresponds to the disease being endemic and is named an endemic equilibrium. To identify the critical threshold for the parameters, we take ×(2.2b)+

provided a key quantity σ ≤ 1, where

is called the modified contact number [23] . The global stability of P 0 when σ ≤ 1 follows from LaSalle's invariance principle. The following theorem (see [23] ), a standard type in mathematical epidemiology, shows that the basic number σ determines the long-term outcome of the epidemic outbreak.

Theorem 2.1. 1. If σ ≤ 1, then P 0 is the unique equilibrium and globally stable in Σ. 2. If σ > 1, then P 0 is unstable and the system is uniformly persistent in Σ 0 .

By this theorem, the disease-free equilibrium P 0 is globally stable in Σ if and only if σ ≤ 1. Its epidemiological implication is that the infected fraction of the population vanishes in time so the disease dies out. In addition, the number σ serves as a sharp threshold parameter; if σ > 1, the disease remains endemic.

In the epidemic literature, another threshold quantity is the basic reproduction number

which is the product of the contact rate β and the average infectious period 1 γ+µ . It is a parameter well known for quantifying the epidemic spread [10, 35] . On the other hand, the contact number σ is defined at all times. In general, we have

Note that for most models, σ = R 0 , both quantities can be used interchangeably. This is the case in our experiments when taking b = 0.

Remark 2.2. Extensive evidence shows that the disease spread rate is sensitive to at least three factors: (1) daily interactions, (2) probability of infection, and (3) duration of illness. The above assertion shows that making efforts to decrease R 0 is essential for controlling propagation of epidemics. Hence measures such as social distancing (selfquarantine, physical separation), washing hands and wearing face coverings, as well as testing / timely-hospitalization can collectively decrease R 0 .

Remark 2.3. The above model if further simplified will reduce to the SIR model or SIS model, which is easier to analyze [13, 5, 34] . It can also be enriched by dividing into different groups [18, 32] , or by considering the spatial movement effects [14, 16] . In this work, we use the SEIR model as a base for our analysis, prediction, and control.

A simple control. In compact form, the SEIR system may be written aṡ

are the column vectors of the state and parameters, respectively. Let g(·) be a loss function (to be specified later) at the final time T , then the problem of determining the model parameters based on this final cost can be cast as a dynamic control problem:

Here the dependence of g on θ is through U (T ). The set of admissible parameters Θ may be estimated from other sources.

For a short time period, the model parameters are near constant, then we can use gradient-based methods to directly update θ, for which ∇ θ g needs to be evaluated. We define the Lagrangian functional

is the Lagrangian multiplier depending on time and can be chosen freely, and U depends on θ through the ODE. A formal calculation gives

Thus ∇ θ g can be determined in the following steps:

• Solve the forward problem for state U ,

• Solve the backward problem for co-state V ,

• Evaluate the gradient of g by

In the context of the optimal control theory [29] , such V exists and is called the co-state function. We note that the (U, V ) dynamical system models the optimal strategies for S, E, I, R populations.

Remark 2.4. It is remarkable that computing gradients of a scalar-valued loss with respect to all parameters is done using only state and co-state functions. without backpropagating through the operations of the solver. In addition, modern ODE solvers may be used to compute both U and V .

2.4. Data-driven optimal control. Now we consider an extended time period with reported data available at intermediate observation times. The data are taken from the reported cumulative infection and death cases. * In such general case the parameters are typically time-varying, we need to derive a more refined data-driven optimal control. Arranging the data in a vector U c = [I c , D c ] at times 0 = t 0 < t 1 < t 2 < ... < t n = T , we aim to * Data used here is publicly available in the CSSEGISandData/COVID-19 GitHub repository, which collects data from official sources and organizations.

• find optimal parameter θ(t) for 0 ≤ t ≤ t n−1 such that the solution to the SEIR system fits the reported data at the grid points {t i } n−1 i=1 as close as possible, and • find desired parameter θ(t) for t n−1 ≤ t ≤ T that is able to control the epidemic spreading at time T at desired values. To achieve these goals, we first define a loss function by

Here λ 1 , λ 2 are user-defined normalization factors. This loss function is composed of two parts: the running cost L which measures the error between the candidate solution (I, D) to the SEIR system and the reported data (I c , D c ) at intermediate times; and the final cost g, which measures the error between the candidate solution (I, D) and the desired data (I d , D d ) at the end time.

We then formulate the task as the following optimal control problem

Motivated by the classical optimal control theory, we derive the necessary conditions for the optimal solution to problem (2.4), stated in the following theorem.

Theorem 2.5. Let θ * be the optimal solution to problem (2.4) and U * be the corresponding state function, then there exists a piece-wise smooth function V * and a Hamiltonian H defined by 

The Pontryagin's maximum principle states the necessary conditions for optimality: assume θ * and U * are the optimal control function and corresponding state of the system, then there exists a function V * and a Hamiltonian H defined for all t ∈ [0, T ] by Notice that the loss function in (2.4) can be rewritten as

where δ(x) is the Dirac delta function, then the corresponding Hamiltonian reads as

Hence (2.8) can be expressed as (2.6), (2.9) has the form of

Integrate (2.11) from t − i − to t + i + and let goes to 0, then for i = n − 1, ..., 1

). Thus in each interval, (2.11) reduces tȯ

which is exactly (2.7). Correspondingly, the Hamiltonian (2.10) reduces to (2.5) for all t ∈ [0, T ] but {t i } n−1 i=1 .

In this section, we present implementation details to solve problem (2.4) via solving the generalized PMP by iteration. We proceed in the following manner. First we make an initial guess θ 0 ∈ Θ. From the control function θ l (t) in the l-th iteration for l = 0, 1, 2, · · · , we obtain θ l+1 (t) in three steps:

Step 1: Solve the forward problem (2.6) to obtain U l .

Step 2: Solve the sequence of backward problems (2.7) to obtain V l .

Step 3: θ l+1 = argmin θ∈Θ H(U l , V l , θ, t) for each t ∈ [0, T ]. This is essentially the method of successive approximations (MSA) [9] . An important feature of MSA is that the Hamiltonian minimization step is decoupled for each t. However, in general MSA tends to diverge, especially if the bad initial guess is taken [9] .

Our goal is to adopt a careful discretization at each step to control possible divergent behavior, instead of simply call an existing ODE solver for Steps 1-2, and an optimization algorithm for Step 3.

Within each interval (t i−1 , t i ), we approximate both U and V at the same m points with step size h = (t i − t i−1 )/m, so that the value of U required in Step 2 are already calculated in Step 1. We use U i,k , V i,k , θ i,k , where i = 1, ..., n and k = 0, ..., m to denote solution values at k-th point in the i-th interval.

Note that population dynamics (births or deaths) may be neglected at a very crude level on the grounds that epidemic dynamics often occur on a faster time scale than host demography, or we can say heuristically that death of an infected individual and subsequent replacement by a susceptible (in the absence of vertical transmission) is equivalent to a recovery event. Hence, in the numerical study, both b and d are taken to be zero.

Below we discuss the discretization of the three iteration steps for the case b = d = 0.

3.1. Forward discretization. We focus on the time interval (t i−1 , t i ). For notation simplicity, we use U k , V k , θ k to represent corresponding values at the k-th point. We discretize the forward problem (2.6) by an explicit-implicit method with the Gauss-Seidel type update:

This gives

where k = 0, 1, ..., m − 1, h = (t i − t i−1 )/m is the step size. The most important property of this update is its unconditional positivity-preserving property, i.e,

irrespective of the size of the step size h. For the starting value in each interval, we have

and notice that

, then (2.7) takes the following forṁ

In a similar fashion to the discretization of U , we discretize this system by an explicitimplicit method by

Then the scheme is given by

/m is the step size. For the starting value in each interval, we have V n,m = ∇ U g(U (T )), V i,m = V i+1,0 + ∇ U L(U (t i )), i = n − 1, · · · , 1. Note that the above discretization is well-defined for any h > 0, hence unconditionally stable since the system is linear in V .

After plugging in U (t), V (t) that is obtained by solving the above forward and backward problems with given θ(t), H may be seen as a functional only with respect to θ. We solve it by the proximal point algorithm (PPA) [30] , that is for any fixed time t ∈ (t i−1 , t i ), i = 1, 2, ..., n,

where l is the index for iteration, τ is the step size. (2) it is numerically stable: PPA has the advantage of being monotonically decreasing, which is guaranteed for any step size τ > 0. In this way, the convergence of minimizing the original loss in (2.4) with a regularization term is ensured, hence address the convergence issue that MSA generally has [1] .

3.4. Algorithms. The above computing process is wrapped up in Algorithm 1. Since we can update {θ i,k } simultaneously, we consider {θ i,k } as a vector and still denote it as θ for notation simplicity. We use · to stand for the L 2 norm.

Require: {U c (t i )} n i=1 : data, t n = T : final time, U 0 : initial data at 0, θ 0 : initial guess, τ : step size for the minimization problem 1: while θ l − θ l−1 / θ l−1 > Tol do 2: for i = 1 to n do 3: for k = 0 to m − 1 do 4:

U i,k+1 ← U i,k (solve the forward problem, refer to (3.1))

U i,0 ← U i−1,m (update the initial condition for ODE solver) 6: V n,m ← ∂ U g(U (T )) (set the initial data for the backward problem) 7: for i = n to 1 do V i,k ← V i,k+1 (solve the backward problem, refer to (3.2)) 10:

) (update the initial condition for ODE solver) 11: θ l+1 ← θ l (solve the minimization problem, refer to (3.5)) 12:

This data-driven optimal control algorithm provides an optimal fitting to both the data and the SEIR model, it can help to reduce the uncertainty in conventional model predictions with standard data fitting.

Initialization: We now discuss how to initialize the control function θ. A simple choice is to take θ i,k 0 to be a constant vector for all 1 ≤ i ≤ n, 1 ≤ k ≤ m, where the value of each component relies on a priori epidemiological and clinical information about the relative parameter magnitude. They vary with the area from where the data were sampled.

We can also use the data {D c (t i )} n i=0 , {I c (t i )} n i=0 to obtain a better initial guess for µ. More precisely, from (2.1), i.e., µ =Ḋ/I we take

The initial data for the forward problem is set as

where N (0) is the initial population of the area analyzed, I c (0), D c (0) are the confirmed infections and deaths on the day of the first confirmed cases, respectively. The initial condition of the backward problem is given in (2.7). In the present setup, the loss function g only depends on I and D, hence the initial condition is given by V n,m = [0, 0, ∂ I g(U (T )), 0, ∂ D g(U (T ))] .

Finally, we should point out that Algorithm 1 with a rough initial guess θ 0 can be rather inefficient when T is large. In such case, we divide the whole interval as 0 = t n 0 < t n 1 < ... < t ns < t n s+1 = T with t n j ∈ {t i } n−1 i=1 for 1 ≤ j ≤ s and apply Algorithm 1 to each subinterval consecutively. This treatment is summarized in Algorithm 2.

Require: {U c (t i )} n i=1 , U 0 : initial data, τ : step size for the minimization problem. Require: {t n j } s+1 j=0 , θ 0 : initial guess of θ on [t n 0 , t n 1 ] 1: for j = 0 to s do 2:

return θ, U Remark 3.3. If parameters vary dramatically, θ(t n j ) may not be a good initialization for θ on [t n j , t n j+1 ]. When this happens we switch to the rough initial guess and test by trial and error. 

We now present experimental results to demonstrate the good performance of our algorithms. In all the experiments, the normalization factors λ 1 , λ 2 in (2.3) are chosen such that the loss with respect to the infection cases is at the same scale as the loss with respect to the death cases. According to the Centers for Disease Control and Prevention (CDC) in the US, the median incubation period is 4-5 days from exposure to symptoms onset, persons with mild to moderate COVID-19 remain infectious no longer than 10 days after symptom onset, therefore, we set

For the step size, considering that β, and γ are almost at the same scale, which is 100 times greater than µ, and the permissible range of β is much large than and γ, we set τ = τ γ = τ , τ β = 100τ and τ µ = τ /100. Figure 1 (a) shows that our data-driven optimal control algorithm learns the data very well. In Figure 1 (b), there is a noticeable peak over the second month, where the value of reproduction number R 0 is very large due to a dramatic increase in infection rate β. After that, R 0 goes down to a lower range, with a slight rise around the 6-th month. Over the last two months we observed another increase in R 0 . Overall, the value of R 0 stays above 1 and the pattern of R 0 is consistent with the increasing trend of the confirmed cases. Short-time prediction. For a short time period, one may expect the transmission to continue the same way; hence the obtained model parameters could be used for prediction over a short coming period. Here, we use the parameters at the 300th day (β = 0.195, = 0.2, γ = 0.1, µ = 1.09 × 10 −4 ) to make prediction for the cumulative number of infection and death cases in the next 14 days. For this prediction, we simply exploit the 4th-order Runge-Kutta solver. Figure 2 shows the prediction result for 300 − 314 days.

Scheduled control. From the above prediction result, we see that without interventions, the amount of confirmed and death cases will increase rapidly. We would like to For instance, starting from the 270-th day, with the goal of controlling the cumulative number of infection and death cases at 9, 746, 063 and 234, 390 on the 300-th day, respectively, we set a pair of values for each day in the 30 days as shown in Figure 3 (a), then learn the parameters from the scheduled data. The results are presented in 3 (b). Figure  3 (a) shows that the goal can be achieved by setting the parameters as what have been learned.

To compare the situations with or without a scheduled control, we also present the reported data and corresponding training results in Figure 3 . In fact, the scheduled intermediate values are obtained by assuming the daily increases were half of the reported daily increases. From Figure 3 (b), we see that the most significant difference occurs in β. This can be roughly interpreted as: if the contact rate β could be reduced by 0.02, the number of confirmed cases over the last 30 days could have been reduced by 50%, though the corresponding R 0 is still greater than 1. For virus propagation to eventually stop, R 0 needs to be less than 1, for which β must be less than γ + µ.

Experimental results for other countries. The coronavirus pandemic continues to affect every region of the world, but some countries are experiencing higher rates of infection, while others appear to have mostly controlled the virus. In order to see the virus dynamics in other regions, we also provide results for some other selected countries Figure 4 and 5, we see that the confirmed cases in the UK and France display similar patterns. Figure 6 shows that China was hit hard early on, but the number of new cases has largely been under control for months.

In this paper, we introduced a data-driven optimal control model for learning the timevarying parameters of the SEIR model, which reveals the virus spreading process among several population groups. Here the state variables represent the population status, such as S, E, I, R, D while the control variables are rate parameters of transmission between population groups. The running cost is of discrete form fitting the reported data of infection and death cases at the observation times. The terminal cost can be used to quantify the desired level of the total infection and death cases at scheduled times. Numerical algorithms are derived to solve the proposed model efficiently. Experimental results show that our approach can effectively fit, predict and control the infected and deceased populations.

The data-driven modeling approach presented in this work is applicable to more advanced models such as with spatial movement effects, interaction of different class of populations, for which we need to formulate mean-field game type models over a spatial domain, see [20] for the mean-field game formulation of a epidemic control problem. On the computational side, our approach involves a non-convex optimization problem, which comes from the multiplicative terms of the SEIR model itself. In future work, we intend to extend our algorithm to more advanced models. 

On the accumulation of perturbations in the linear systems with two coordinates

Population biology of infectious diseases: Part I

Optimal control of deterministic epidemics

Propagation of epidemics along lines with fast diffusion

Mathematical models in population biology and epidemiology

Introduction to the mathematical theory of control

Cancer immunotherapy, mathematical modeling and optimal control

Neural ordinary differential equations, Conference on Neural Information Processing Systems (NIPS)

Method of successive approximations for solution of optimal control problems

Mathematical epidemiology of infectious diseases

Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity

Modeling epidemics with variable contact rates

The mathematics of infectious diseases

Traveling waves for a simple diffusive epidemic model

Optimal control problem of an SIR reaction-diffusion model with inequality constraints

Mathematical models of the spread of infection

A contribution to the mathematical theory of epidemics

Global properties of SIR and SEIR epidemic models with multiple parallel infectious stages

Dynamics and optimal control of a nonlinear epidemic model with relapse and cure

Controlling propagation of epidemics via mean-field games

Optimal control applied to biological models

Dynamic stability of an SIQS epidemic network and its optimal control

Global dynamics of a SEIR model with varying total population size

Global stability for the SEIR model in epidemiology

Dynamical behavior of epidemiological models with nonlinear incidence rates

First-principles machine learning modelling of COVID-19

Dynamic models of infectious diseases as regulators of population sizes

On the optimal control of a deterministic epidemic

The mathematical theory of optimal processes

Monotone operators and the proximal point algorithm

Optimal control of some simple deterministic epidemic models

Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence

Epidemic and demographic interaction in the spread of potentially fatal diseases in growing populations

Exact analytical solutions of the susceptible-infectedrecovered (SIR) epidemic model and of the SIR model with equal death and birth rates

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Given discretized U α and V α with α = {i, k}, we solve (3.3) on grid points, that is to solveSince H is smooth, the above formulation when the constraint is not imposed is equivalent to the following. The special form of H allows us to obtain a closed form solution:. Now taking the constraint θ ∈ Θ into consideration, we simply project the solution (3.5) back into the feasible region after each iteration, that is θ α l+1 = clip(θ α l+1 , Θ, Θ) where Θ, Θ are element-wise lower bound and upper bound of Θ, respectively. That guarantees the output is constrained to be in Θ.