Robust optimal control of compartmental models in epidemiology: Application to the COVID-19 pandemic.

In this paper, a spectral approach is used to formulate and solve robust optimal control problems for compartmental epidemic models, allowing the uncertainty propagation through the optimal control model to be represented by a polynomial expansion of its stochastic state variables. More specifically, a statistical moment-based polynomial chaos expansion is employed. The spectral expansion of the stochastic state variables allows the computation of their main statistics to be carried out, resulting in a compact and efficient representation of the variability of the optimal control model with respect to its random parameters. The proposed robust formulation provides the designers of the optimal control strategy of the epidemic model the capability to increase the predictability of the results by simply adding upper bounds on the variability of the state variables. Moreover, this approach yields a way to efficiently estimate the probability distributions of the stochastic state variables and conduct a global sensitivity analysis. To show the practical implementation of the proposed approach, a mathematical model of COVID-19 transmission is considered. The numerical results show that the spectral approach proposed to formulate and solve robust optimal control problems for compartmental epidemic models provides healthcare systems with a valuable tool to mitigate and control the impact of infectious diseases.

Introduction

Optimal control is a research field with many applications in economics [1], logistics [2], transportation [3], engineering [4], biology [5], and biomedicine [6]. It deals with the optimization of some performance index imposed on an underlying dynamical system subject to constraints. In this framework, controls are functions of time that describe allowable external influences on the system that induce a system response and, based on this response, an objective functional is evaluated that is considered as a measure of the performance of the system’s behavior. Optimal control addresses the question of finding the controls that steer the dynamical system from an initial to a final state while optimizing this objective functional.

External disturbances and model uncertainties affect the solutions of optimal control problems. However, they are difficult to include in a deterministic mathematical description of the system. Stochastic models must be used, in which the dynamical system is described by a set of stochastic Ordinary Differential Equations (ODE), which, in this paper, are represented by a set of ODE with random parameters. Optimal Control Problems (OCP) for this class of systems, called stochastic dynamical systems, are usually referred to as stochastic optimal control problems. Stochastic optimal control problems require robust solutions to be found, that is, solutions in which the control task is accomplished independently of the realization of the random parameters. This entails including statistical moments and probabilities of events in the constraints or the objective functional of stochastic optimal control problems, converting them into robust optimal control problems.

There are different methods for dealing with robust OCP, such as integral averaging [7], adjustable uncertainty sets [8], adaptive dynamic programming [9], and neural networks [10]. In this paper, a spectral approach is employed to formulate and solve robust OCP for compartmental epidemic models, which allows the uncertainty propagation through the optimal control model to be represented by means of a polynomial expansion of its stochastic state variables. In particular, a Polynomial Chaos Expansion (PCE) [11] is used to convert the robust OCP with stochastic ODE and stochastic constraints into an equivalent deterministic OCP with deterministic ODE and deterministic constraints. The PCE is a probabilistic method that consists in the projection of the optimal control model on a basis of orthogonal stochastic polynomials in the random parameters. This stochastic projection provides a compact and efficient representation of the variability of the optimal control model with respect to the random parameters.

The mitigation and control of the spread of an epidemic disease has all the characteristics of an OCP. The aim is to minimize some objective functional usually related to the number of infectious individuals or the cost of the pharmacological control measures of the disease, while the underlying dynamical system describes the progression of the disease transmission and its interactions with the mitigation and control strategies. A paradigmatic case is to find the optimal vaccination schedule when vaccination supplies are limited. In this paper, an objective functional that combines the number of infected individuals with the vaccination cost is considered in the presence of an upper bound on the daily number of vaccinated people expressed as a percentage of the susceptible individuals. The goal is to determine the optimal vaccination strategy.

One of the main approaches to the mathematical modeling of epidemic diseases consist in describing the disease transmission by means of compartmental models, in which the population is subdivided into compartments [12]. In the last two decades, the analysis of the evolution of different epidemic outbreaks has relied on these models [13], [14], [15], [16]. More recently, various research works based on compartmental epidemic models have been carried out to study different aspects of the spread, mitigation, and control of the COVID-19 disease. They include a classification of models [17], an estimation of the size of the epidemic [18], a study of the effects of pharmacological measures [19], and the analysis of the spread of the COVID-19 disease focused on the decrease in susceptible individuals [20], the impact of superspreaders  [21], and the different non-pharmaceutical containment and mitigation strategies, such as testing, contact tracing, and social distancing [22], [23].

The optimal control approach has been widely used to study compartmental epidemic models. For the sake of brevity, the following recent works are highlighted. A nonlinear epidemic model on complex networks is introduced in [24] in order to better understand and utilize the quarantine control when encountering outbreaks of infectious diseases. In [25], an isolation problem for a specific epidemic model describing the spread of AIDS is explicitly solved by means of optimal impulse control. A reaction–diffusion model is studied in [26] to derive an efficient vaccination strategy for influenza outbreaks. In [27], a COVID-19 disease spread model with free terminal optimal time is defined in which the goal is to reduce the size of susceptible, infected, exposed, and asymptomatic groups in order to consequently eradicate the infection by using quarantine and treatment of infected people. In [28], the effects of non-pharmacological strategies such as quarantine, isolation, and public health education are studied as time-dependent interventions to ascertain their contributions to the transmission dynamics of COVID-19. Optimal control strategies for the Zika infection with mutations that cause defects in the newly birth for infected pregnant woman are formulated in [29].

Fewer works have been devoted to the robust optimal control of compartmental epidemic models. In [30], the evolutionary dynamics of HIV in the presence of drug therapy is modeled by means of a class of positive systems in which the control signal enters bilinearly with the state. In [31], a general heterogeneous susceptible–infected–susceptible model is considered by taking into account noisy transition rates. In [32], the optimal policy for controlling the COVID-19 epidemic dynamics using both lockdown and detection intervention levers is identified, which takes into account the trade-off between the sanitary and the socio-economic cost of the pandemic, together with the limited capacity of intensive care units. In [33], a co-infection model is developed for human papillomavirus and Chlamydia trachomatis with cost-effectiveness optimal control analysis.

In this paper, the problem of finding optimal control strategies for compartmental models in epidemiology is studied. The problem is modeled as a robust optimal control problem, and a spectral approach is used to represent the uncertainty propagation through the optimal control epidemic model using its PCE. More specifically, the statistical moment-based PCE introduced in [34] is employed as a surrogate model of the compartments of the epidemic model, in which a spectral expansion of the stochastic state variables is used, which permits the computation of their main statistics to be carried out in terms of the deterministic state variables. The availability of these statistics allows the robust formulation of the stochastic OCP to be solved. This problem is reformulated as a deterministic OCP with deterministic ODE and deterministic constraints, and is numerically solved using a direct transcription technique. More specifically, the Hermite–Simpson collocation method [35] is employed. The resulting Nonlinear Programming (NLP) problem is solved using the NLP solver IPOPT [36] with the AMPL interface.

The proposed robust formulation of the OCP provides the designers of the strategy of the optimal control epidemic model the capability to increase the predictability of the results by simply adding upper bounds on the variability of the state variables. Furthermore, the probability distributions of the stochastic state variables can be efficiently estimated and a global sensitivity analysis, based on the variance of the random parameters of the optimal control model, can be conducted, with the goal of determining which random parameters have the greatest impact on the variability of the state variables at each time instant and at a low computational cost.

To show the effectiveness of the proposed approach, the mathematical model of the SARS-CoV-2 virus transmission dynamics introduced in [37] and modified in [38] is considered, which includes a vaccination strategy. In particular, assuming the implementation of a vaccination plan with a constraint on the daily available vaccines, the effects of the uncertainty in the implementation of social distancing measures and the testing of susceptible individuals are quantified with the aim of finding the optimal vaccination strategy. The results of various numerical experiments show that the spectral approach to the robust OCP formulation of compartmental epidemic models provides healthcare systems with a valuable tool to mitigate and control the impact of infectious diseases.

This paper is organized as follows. In Section 2, starting from the definition of a continuous deterministic OCP, in which the presence of uncertainties in some of its parameters is assumed, the general formulation of the robust OCP is presented. The PCE technique, which is used to represent the propagation of these uncertainties through the optimal control model, is then introduced in Section 3. The methodology employed to compute the main statistics, reformulate the robust OCP, estimate the distributions of the stochastic state variables, and conduct a global sensitivity analysis is described in Section 4. The reformulation of the robust OCP as a deterministic OCP through the PCE technique is presented in Section 5. The mathematical model of the SARS-CoV-2 virus transmission, used to show the effectiveness of the proposed procedure, is described in Section 6. The robust optimal control of COVID-19 transmission dynamics with mixed state-control constraint is formulated in Section 7. In Section 8, the results of several numerical experiments are reported and analyzed. Finally, some conclusions are drawn in Section 9.

Robust optimal control

In this section, following [39], the formulation of the robust OCP considered in this paper is presented, which is based on the definition of a deterministic continuous OCP.

Statement of the deterministic optimal control problem

The deterministic OCP considered in this paper can be stated as follows: where denotes time, is the vector of state variables, and is the vector of control variables. The initial and final times are denoted by and , and and represent the vector of initial and final states, respectively. The objective functional in (1a) is given in Bolza form, that is, it is expressed as a combination of a Mayer term , which represents a terminal cost, and a Lagrange term , which represents a running cost. The set of differential equations that represents the dynamical system is described by Eq. (1b), Eq. (1c) denotes the so-called mixed state-control constraints, and the initial and final conditions, that is, the boundary conditions, are represented by Eq. (1d).

Note that, without loss of generality, the objective functional can be formulated in Mayer form, since the Lagrange term can be transformed into a Mayer term by introducing a new state variable, , and adding a new differential equation together with its corresponding initial condition, as follows: Then, the Lagrange term in the objective functional (1a) can be replaced by the term .

Numerical resolution of the deterministic optimal control problem

In general, it is very difficult to find an analytical solution to a nonlinear OCP. Thus, the common practice is to use numerical methods to obtain solutions. In this work, a direct numerical method is employed to transcribe the deterministic OCP defined in (1) into a NLP problem. More specifically, the Hermite–Simpson direct collocation method is used [35]. The set of constraints of the resulting NLP problem includes the Hermite–Simpson system constraints that correspond to the differential constraint (1b) and the discretized versions of the other constraints of the OCP. They include the mixed state-control constraints (1c) and the boundary conditions (1d). To solve the resulting NLP problem, the open source IPOPT solver [36] with the AMPL interface has been employed.

Statement of the stochastic optimal control problem

Let be a probability space, where is the space of events, is a -algebra, and is a probability measure. Then, if uncertainties are to be taken into account, a probabilistic version of the deterministic OCP defined in (1) must be considered, which can be stated as follows: where denotes the vector of random parameters that describe the uncertainties, and the rest of the functions and variables are defined as in (1). Thus, the objective functional, the differential equations, and the constraints functions in (3) are stochastic functions. Note that, in this formulation, the control variables remain deterministic, while the state variables are functions of both the time and the vector of random variables , and the stochastic relations are assumed to be satisfied almost surely (a.s.).

Statement of the robust optimal control problem

A robust version of the OCP stated in (1) can be derived from the stochastic OCP defined in (3). On the one hand, the robustness of the objective functional (1a) can be achieved by considering a weighted sum of the mean value and the standard deviation of the stochastic objective functional (3a). On the other hand, the robustness of the mixed state-control constraints (1c) and the boundary values (1d) can be modeled by splitting the mean and standard deviation conditions. Thus, the mean values of the stochastic constraints (3c) and the boundary conditions (3d) are forced to satisfy the same requirements stated in the deterministic OCP defined in (1), whereas the standard deviation of these stochastic constraints and boundary conditions are forced to satisfy the upper bounds specified by the designers of the optimal control model.

Thus, the robust OCP considered in this work can be stated as follows: where and denote the expected value and the standard deviation operators, respectively, and the remaining variables and functions are defined as in (1) and (3). As mentioned above, the constant coefficient in (4a) and the constant upper bounds , , and for the standard deviations in constraints (4c) and boundary conditions (4d), must be specified by the designers of the optimal control model.

Moment-based arbitrary polynomial chaos expansion

In this section, following [34], the technique used to model the propagation of the uncertainties through the stochastic OCP defined in (3) is presented.

Polynomial chaos expansion

Let be a vector of independent random variables in the probability space introduced in Section 2.3. Then, a multi-dimensional polynomial expansion dependent on the random vector can be defined to represent the vector of stochastic state variables of the stochastic OCP defined in (3).

In particular, from a theoretical point of view, the vector can be expressed as a linear combination of stochastic multivariate orthonormal polynomials , such that with deterministic coefficients  [11].

In practice, from a computational point of view, a truncated version of the expansion (5) must be considered. Thus, vector can be approximated by a linear combination of stochastic multivariate orthonormal polynomials as follows: Each of the multivariate orthonormal polynomials , , in (6) can be built as a product of univariate orthonormal polynomials , , , where index refers to the random variable and index refers to the order of the univariate orthonormal polynomial.

For , each polynomial of the expansion (6) can be calculated as where the index matrix represents a graded lexicographic ordering, that is, the rows of the index matrix indicate which order of each univariate polynomial contributes to a particular multivariate polynomial. Note that, if a full tensor product is used, the number of terms in (6), for an expansion of order , is

The orthogonality condition implies that, for each univariate random variable , , the corresponding univariate orthonormal polynomials for with must fulfill with the Kronecker delta and where the th orthonormal polynomial for each random variable is defined by means of its coefficients as

Due to the orthonormality of the polynomials, the coefficients , , of the expansion (6) can be obtained as In practice, integral (7) can be solved by means of Galerkin projection, collocation, or numerical integration, among others. In this paper, a numerical integration based on the Gaussian quadrature rule has been employed.

Gaussian quadrature rule based on statistical moments

A Gaussian quadrature rule is a well-known approach that provides the collocation points and weights , , for the numerical integration of a univariate function defined on a measure space and, thus, on a probability space . This Gaussian numerical integration gives an approximation of an integral of the form Note that, for the sake of clarity of exposition, in this section the dependency of function on time has been omitted.

In case of multiple random variables, multivariate quadrature formulas can be obtained from one dimensional quadrature rules. The most natural and straightforward extension of the one dimensional quadrature rules to the multivariate quadrature case is the full tensor product quadrature rule. Given a vector of independent random variables , the full tensor product quadrature formula for a multivariate function can be expressed as

Note that, on the one hand, the optimal values of the collocation points and weights, and their corresponding orthogonal polynomials basis, can be obtained for a large number of parametric distributions using the generalized Polynomial Chaos (gPC) expansion method introduced in [11]. On the other hand, a more general approach is provided by the arbitrary Polynomial Chaos (aPC) expansion [40], which allows different kind of descriptors to be considered, including data sets without a known parametric probability distribution.

Given a random variable , in the methodology proposed in [40], the PCE is obtained from the Hankel matrix of moments, which can be defined as where the entry , , denotes the th raw statistical moment of random variable . Since it is a positive definite matrix, a Cholesky decomposition of the Hankel matrix of the form can be computed, where Let the inverse of matrix be Then, as shown in [41], the entries of matrix form an orthogonal system of polynomials such that

In order to avoid the inversion of the matrix , explicit analytic formulas can be derived to obtain the polynomial coefficients of the orthogonal polynomials directly from the entries of the matrix . In particular, the elements can be used to determine the coefficients and of the following three-term recurrence relation where and . More specifically, coefficients and can be expressed in terms of entries as follows with and .

The optimal collocation points and weights for the orthogonal polynomials can be efficiently calculated using the coefficients and . The above three-term recurrence relation can be rewritten by means of the positive definite symmetric tri-diagonal Jacobi matrix In this way, on the one hand, the roots of the polynomial of order are obtained directly as the eigenvalues , , of matrix . On the other hand, the weights can be calculated as where , , represents the first component of the normalized eigenvector corresponding to the th eigenvalue of matrix .

Note that, when the number of random variables grows, the full tensor Gaussian quadrature rule (8) becomes expensive from a computational point of view. The usual approach to overcome this drawback is to consider quadrature rules based on sparse grids, which reduce significantly the number of nodes and weights, while keeping a reasonable accuracy level. The most common sparse grid quadrature rule is the Smolyak quadrature rule based on sparse tensor product spaces [42].

Uncertainty quantification

One of the main advantages of the PCE approach is its computational efficiency in both calculating the statistics and distributions of the optimal control model (3) and analyzing its sensitivity to the presence of uncertainties through the surrogate model (6).

Statistics of the state variables

Given the multi-dimensional polynomial expansion (6), the statistics of the stochastic vector of state variables can be expressed in terms of the coefficients in a straightforward manner. For example, the mean and variance of the stochastic vector of state variables can be calculated using the following expressions

Moreover, the quadrature rule described in Section 3.2 provides another simple and computationally efficient way to calculate the statistics of by means of scalar products of vectors. In particular, the statistics of based on the computation of its moments can be efficiently computed through the integration rule (8) since, for every , the following equation holds For instance, the mean and variance of the stochastic vector of state variables can be obtained as with , and denoting the vector of Gaussian quadrature weights.

Distributions of the state variables

The surrogate model provided by the multi-dimensional polynomial expansion (6) can be employed to estimate, at a low computational cost, the PDF of the stochastic vector of state variables . First, a data smoothing problem is formulated and then a kernel density estimation approach [43] is conducted, in which the data are extracted from the surrogate model. The Cumulative Density Function (CDF) of the stochastic vector of state variables can be estimated in a similar way.

Sensitivity analysis

The multi-dimensional polynomial expansion (6) can be also employed to perform a global sensitivity analysis based on the variance of the stochastic vector of state variables , with the aim of quantifying how much of the variance of is explained by the variance of each random variable. Under the assumptions that the stochastic vector of state variables is square integrable and the components of the vector of random variables are independent, this analysis is based on the calculation of the so-called Sobol indices [44]. Sobol’ indices can be also directly obtained from the coefficients of the aPC expansion (6). More specifically, the Sobol’ index for the th component of the vector of random variables can be computed, following [45], as where is the variance calculated in (9) and denotes the set of multi-indices such that the computation of the Sobol’ index only includes terms that depend on the random variable , , namely,

Reformulation of the robust optimal control problem

Following [39], the numerical resolution of the robust OCP stated in (4) is based on the reformulation of the objective functional (4a), the stochastic differential Eqs. (4b), the constraints (4c), and the boundary conditions (4d) in terms of the aPC expansion described in Section 3.

Let , , be the set of nodes provided by the aPC expansion of order of the stochastic state vector with respect to the random vector . Then, the robust OCP defined in (4) can be rewritten as follows: where , , and , , are the coefficients of the aPC expansion of the state variable introduced in (2).

Note that from (9), (10) it is easy to see that the terms , , , and only depend on the deterministic state vector . Thus, the optimal control problem (11) is a deterministic OCP, which can be numerically solved by using the techniques described in Section 2.2. More specifically, the computational strategy employed to calculate the optimal controls and state variables of the robust OCP (4) is schematized in Fig. 1.

Fig. 1

Computational strategy for the resolution of the robust OCP stated in (4).

Mathematical model of COVID-19 transmission

The availability of various effective COVID-19 vaccines, which represent an effective pharmacological measure against this pandemic, permits new strategies to be implemented, in which non-pharmaceutical public health measures and vaccination are combined. A mathematical model of the SARS-CoV-2 virus transmission dynamics with vaccination strategy is considered in this paper, which takes into account the testing of susceptible individuals and the implementation of social distancing measures. There are inherent uncertainties in the sensitivity and specificity of COVID-19 tests, as well as about the degree of application of social distancing measures, that should be taken into consideration in order to obtain a more realistic model. However, the presence of uncertainties in the application of these non-pharmacological strategies results in a complex stochastic nonlinear model that is difficult to solve.

A modified SEIR model, whose parameters are adjusted using data from the World Health Organization (WHO), is considered in this work following [37]. This model, schematically shown in Fig. 2, is denoted as SEIsIaQR and has 6 compartments. It has been shown in [37] that the SEIsIaQR model is non-chaotic and asymptotically stable and has a unique equilibrium state. It is represented by the following set of Ordinary Differential Equations (ODE): Note that, following [38] a vaccination term has been included as a control variable. Thus, the dynamical system (12) has 7 state variables, 1 control input, and 8 parameters.

Fig. 2

Schematic representation of the SEIsIaQR compartmental model with vaccination strategy [38].

The state variables of the dynamical system are , , , , , and , which represent the number of susceptible, exposed, symptomatic, asymptomatic, isolated, and recovered individuals, respectively, and , which represents the total population. Thus, the following condition must be satisfied: with denoting the final time considered in the ODE system.

The control variable of the system is , which represents the daily vaccination rate. It ranges from 0 to , where coefficient is an upper bound that must be specified by the designers of the optimal control model. Thus, the vaccination strategy in model (12) is defined by the term .

The parameters in model (12) that are independent of the mitigation strategies are the mean serial time interval, , the mean incubation period, , the infectious period, , the ratio between infectiousness of asymptomatic and symptomatic individuals, , and the population ratio that remains asymptomatic, . These parameters have been set at the same values as in [38].

The parameters in model (12) related to the use of non-pharmacological mitigation measures are the replication factor, , the isolation rate of symptomatic individuals, , and the identification and isolation rate of asymptomatic individuals, . Thus, assuming the availability of an effective vaccine and following [38], in this paper, an immunization term is added to the model proposed in [37], a robust OCP is formulated, and the effects on the disease transmission of the uncertainty in the mitigation measures represented by parameters , , and are quantified. In addition, surrogate models for the compartments of the model (12) are computed, which are used to estimate the probability distributions of the stochastic state variables. Furthermore, a global sensitivity analysis is performed to determine the random parameters that have the greatest impact on the variability of the stochastic state variables.

Robust optimal control of COVID-19 transmission

Depending on the kind of objective functional considered and the constraints imposed, different strategies to control the spread of COVID-19 represented by the model (12) can be defined. According to [46], the objective functional to be minimized has been assumed to be a combination of the number of infectious individuals and the overall cost of the vaccination during a fixed time period, that is, where is a weighting parameter that must be specified by the designers of the optimal control model. Note that, larger values of give less importance to the vaccination cost. Additionally, following [47], an upper bound on the daily available vaccines has been imposed. This limitation on the number of daily available vaccines must be introduced when the vaccination time interval is large enough or the capability to vaccinate each day may dictate the need for a constraint on the number of daily vaccinated people. Thus, the following constraint has been considered where parameter denotes an upper bound on the number of daily vaccinated people. Note that expression (14) is a mixed state-control constraint of the form (1c). Therefore, the robust optimal control of the SEIsIaQR compartmental model (12) with mixed state-control constraint can be stated as follows: where is defined as in (13), , and .

Numerical experiments

Several numerical experiments are conducted to illustrate the application of the proposed spectral approach to the resolution and uncertainty quantification of the robust OCP defined in (15). According to the numerical experiments presented in [38], a non-dimensionalized version of the system (12) is considered in all these experiments, namely, , and the following initial conditions are assumed at time instant : Moreover, it is assumed that, in order to develop immunity, two doses of vaccine are required and 95% vaccine efficacy, with an upper bound of on the proportion of susceptible daily vaccinated people.

Case Study 1

In this case study, the upper bound on the daily vaccination rate is set to and the numerical experiments are solved assuming that social distancing actions are barely applied and the number of daily tests allow both infected asymptomatic and symptomatic individuals to be identified and isolated. In particular, according to [38], is assumed to follow a gamma distribution G(3500, 0.001), a beta distribution B(160, 10), and a beta distribution B(160, 160), as shown in Fig. 3. Then, these three random parameters are modeled using an aPC expansion of order , the constant bounds for the standard deviations in the boundary conditions (15k)–(15l) are set to , and the weighting parameters and are set to 1 and 0.01, respectively.

Fig. 3

PDF of the random parameters related to the use of non-pharmacological mitigation strategies.

Three solutions for the robust OCP (15) have been computed using three different values for the upper bound . The mean value of the state variables that describe each compartment of the model obtained in these solutions are represented, together with a 95% confidence envelope, in Fig. 4, whereas the corresponding control variables, which represent daily vaccination rates, are depicted in Fig. 5. The curves for , , and are shown in red, blue, and green, respectively. For the sake of comparison, Fig. 4, Fig. 5 also show in black dashed lines the solution obtained by solving the deterministic version of the robust OCP (15), in which the random parameters have been replaced by their mean values.

Fig. 4

Case Study 1. Mean values of the proportions of individuals in each compartment obtained for different values of the upper bound , together with the corresponding 95% confidence envelopes. The curves for , , and are shown in red, blue, and green, respectively. Proportions in dashed black lines correspond to the solution of the associated deterministic OCP.

Fig. 5

Case Study 1. Daily vaccination rate obtained for different values of the upper bound . The curves for , , and are represented in red, blue, and green colors, respectively. The vaccination rate in dashed black line corresponds to the solution of the associated deterministic OCP.

Fig. 4 shows that there are significant differences in the solutions obtained using the three different values for . These differences are noticeable in both the shape of the curves and the dispersion of their envelopes. Lower values of correspond to lower peaks and less dispersion of the corresponding curves. For example, Table 1 shows the 95% confidence intervals for the maximum percentages of individuals in compartments , , , and . It is worth noting that, in the worst-case scenario, the percentage of symptomatic individuals in a single day rises to 0.0290%, 0.0144%, and 0.0090% for , , and , respectively. While the percentage in the first case could place considerable stress on the healthcare system, the last two cases could be managed by a sufficiently resilient healthcare system.

Table 1

95% confidence intervals for the maximum percentages of individuals in compartments , , , and that correspond to different values of the upper bound .

	E	I_s	I_a	Q
ɛ_S=10^−3	(0.3006, 1.0356)%	(0.0085, 0.0290)%	(0.0445, 0.1777)%	(0.2618, 0.8037)%
ɛ_S=10^−4	(0.2166, 0.5141)%	(0.0059, 0.0144)%	(0.0324, 0.0882)%	(0.1795, 0.3963)%
ɛ_S=10^−5	(0.1970, 0.3197)%	(0.0056, 0.0090)%	(0.0301, 0.0549)%	(0.1564, 0.2407)%

Furthermore, this solution implies that (16.93, 17.59)%, (31.74, 32.55)%, and (41.42, 41.83)% of the susceptible individuals become immune through vaccination and, after six months, the percentage of immunized individuals is , , and for , , and , respectively. This last estimation would be close to, but fall just short of, 66.7%, the estimated threshold for achieving herd immunity [48].

It is worth mentioning that the solution obtained by solving the deterministic version of the robust OCP (15) can be regarded as a reference solution, which provides a tool to analyze the effects of the upper bound on the behavior of the optimal control and state variables. It can be seen in Fig. 4, Fig. 5 that, assuming , the optimal control and the mean of the optimal state variables obtained by solving the robust OCP are very similar to those obtained by solving the deterministic OCP. Thus, the robust constraint on the variance of the sate variable in (15h) is irrelevant if the value of the upper bound is greater than or equal to 10−3.

In Fig. 6, the solution obtained using the proposed robust OCP-based methodology is compared to other two solutions. First, the ODE system (12) has been solved assuming a constant vaccination rate , for all , and supposing that the parameters , and are the same random parameters as those considered in the statement of the robust OCP. The corresponding results, which are non-optimal solutions, are represented in green lines. Then, the same stochastic ODE system has been solved removing the vaccination term , namely setting the control input to zero. The corresponding state variables of the uncontrolled system are represented in red lines. The optimal state variables obtained by solving the robust OCP with are represented in blue lines.

Fig. 6

Case Study 1. Mean values of the proportions of individuals in each compartment, together with the corresponding 95% confidence envelopes, obtained by solving the robust OCP with (blue), the ODE system with constant control (green), and the uncontrolled ODE system (red).

It can be seen in Fig. 6 that the implementation of a non-optimal vaccination strategy would imply a significant increase in both the proportion of symptomatic and asymptomatic infected individuals and the dispersion of their corresponding curves. Moreover, in the absence of a vaccination strategy, a noticeable outbreak of the infectious disease would occur, which could place considerable stress on the healthcare system.

Consider variable shown in Fig. 7, obtained in the solution to the numerical experiment with . This variable represents the proportion of symptomatic individuals provided by the solution to this experiment. A Gaussian kernel estimator is employed to approximate the distributions of at time instants days, days, and days. Fig. 8 shows the estimated PDF and CDF of these three random variables. Many different statistical analyses can be conducted using these distributions. For example, at time instant days, the probability that the percentage of symptomatic individuals will exceed 0.015% of the population, is 0.0200, whereas this probability rises to 0.7373 if the experiment is solved assuming .

Fig. 7

Case Study 1. Time slices days (red), days (orange), and days (green) of variable , and the corresponding 95% confidence envelope in the solution of the experiment with .

Fig. 8

Case Study 1. PDF and CDF of variable at time instants days (red), days (orange), and days (green) in the solution to the experiment with .

Similarly, a multivariate statistical analysis can be carried out. For instance, the correlation matrix for the distributions of symptomatic individuals at time instants days, days, and days, is Thus, despite their separation in time, these three distributions show a strong linear association. Moreover, the estimated joint PDF and CDF of variable at time instants days and days, which are shown in Fig. 9, allow conditional probabilities to be calculated. For instance, assuming that at time instant days, the observed percentage of symptomatic individuals of the population exceeds 0.008%, the probability that at time instant days, the percentage of symptomatic individuals exceeds 0.015% would be 0.9073.

Fig. 9

Case Study 1. Joint PDF and CDF of variable obtained in the solution to the experiment with at time instants days and days.

The Sobol’ indices of all the state variables obtained in the solution to the experiments are plotted in Fig. 10, where the indices computed assuming that , , and are shown using dashed, solid, and dotted lines, respectively. The Sobol’ indices corresponding to the random parameters , , and are plotted in blue, orange, and green colors, respectively. Note that, for all the state variables, the proportion of variance explained by the variance of parameter is negligible, whereas between 60 and 80 percent of the variance is due to the variance of parameter . Moreover, for all the state variables, there are no significant differences between the Sobol’ indices computed from the models using the three different values of the upper bound . Thus, the proportion of variance of the state variables due to the variance of each random parameter does not change when the upper bound changes.

Fig. 10

Case Study 1. Sobol’ indices for the numerical experiments with (dashed), (solid), and (dotted). The corresponding Sobol’ indices for the random parameters , , and are plotted in blue, orange, and green, respectively.

Case Study 2

In this case study, the effects of a decrease in the rate of identification and isolation of asymptomatic individuals, , on the behavior of the optimal control and state variables is analyzed. In particular, is assumed to follow a beta distribution B(40, 160) instead of a beta distribution B(160, 160), as shown in Fig. 11. Moreover, the upper bound on the daily vaccination rate is set to . The rest of random and deterministic parameters are assumed to take the same values as in Case Study 1.

Fig. 11

PDFs of the identification and isolation rate of asymptomatic individuals, , for Case Study 1 (right) and Case Study 2 (left).

Three solutions for the robust OCP (15) have been computed using three different values for the upper bound . The mean value of the state variables that describe each compartment of the model obtained in these solutions are represented, together with a 95% confidence envelope, in Fig. 12, whereas the corresponding control variables, which represent daily vaccination rates, are depicted in Fig. 13. The curves for , , and are shown in red, blue, and green, respectively.

Fig. 12

Case Study 2. Mean values of the proportions of individuals in each compartment obtained for different values of the upper bound , together with the corresponding 95% confidence envelopes. The curves for , , and are shown in red, blue, and green, respectively.

Fig. 13

Case Study 2. Daily vaccination rate obtained for different values of the upper bound . The curves for , , and are represented in red, blue, and green colors, respectively.

It can be seen in Fig. 12 that a huge outbreak of the infectious disease occurs when the upper bound in (15h) is set to , which would imply an overwhelmed healthcare system and a dramatic proportion of deaths. This is due to the fact that the corresponding optimal vaccination strategy shown in Fig. 13 is not enough to make up for the decrease in the rate of identification and isolation of asymptomatic individuals. Therefore, in order to achieve a satisfactory optimal vaccination strategy, it is necessary to find a small enough threshold for the dispersion of the state variable . Indeed, it can be seen in Fig. 12 that the optimal vaccination strategies obtained for and , which grow at a high rate during the first 60–80 days, are able to compensate the decrease in the rate of identification and isolation of asymptomatic individuals. However, to obtain these optimal vaccination profiles it is necessary to increase the upper bound on the daily vaccination rate . In fact, the robust OCP (15) is infeasible if is assumed to follow a beta distribution B(40, 160) but the upper bound on the daily vaccination rate remains the same as in Case Study 1, namely . These results supports the hypothesis that the early identification and isolation of asymptomatic individuals may reduce the undesirable effects of the pandemic, while saving the amount of administered vaccines.

Case Study 3

In this case study, the sensitivity of the robust OCP (15) to the weighting parameter of the objective functional (13) is analyzed. More specifically, three solutions of the robust OCP described in the Case Study 1 are computed using three different values of the weighting parameter , namely 0.01, 0.1, and , and assuming that the upper bound in (15h) is set to .

The mean value of the state variables that describe each compartment of the model obtained in these solutions are represented, together with a 95% confidence envelope, in Fig. 14, whereas the corresponding control variables, which represent daily vaccination rates, are depicted in Fig. 15. The curves for , , and are shown in blue, red, and green, respectively.

Fig. 14

Case Study 3. Mean values of the proportions of individuals in each compartment obtained for different values of the weighting parameter , together with the corresponding 95% confidence envelopes. The curves for , , and are shown in blue, red, and green, respectively.

Fig. 15

Case Study 3. Daily vaccination rate obtained for different values of the weighting parameter . The curves for , , and are represented in blue, red, and green colors, respectively.

Fig. 14 shows that there are significant differences in the solutions obtained using the three different values of the weighting parameter . As already mentioned in Section 7, larger values of give less importance to the vaccination term and more importance to the number of infectious individuals in minimizing the objective functional (13). Indeed, it can be seen in Fig. 14 that larger values of correspond to lower peaks of the curves of both infected symptomatic and asymptomatic individuals. Moreover, the larger the value of , the lower the dispersion of the envelopes of the corresponding curves.

Computational aspects

The numerical experiments have been carried out on a desktop computer with a 4 GHz i7 processor with 16 GB RAM under Mac OS X 10.12 operating system. Each one of these numerical experiments consists of 39241 variables, 38880 equality constraints, and 540 inequality constraints.

Table 2 shows, for each numerical experiment, the computation time along with the number of iterations (#ite), the number of objective function evaluations (#ofe), the number of objective gradient evaluations (#oge), the number of equality constraint evaluations (#ece), the number of inequality constraint evaluations (#ice), the number of equality constraint Jacobian evaluations (#ecje), the number of inequality constraint Jacobian evaluations (#icje), and the number of Lagrangian Hessian evaluations (#lhe).

Table 2

Summary of the computational time and the iteration and function evaluation counts for each numerical experiment.

		ct	#ite	#ofe	#oge	#ece	#ice	#ecje	#icje	#lhe
Case Study 1	ɛ_S=10^−3	102.16 s	35	48	36	48	48	36	36	35
	ɛ_S=10^−4	106.51 s	36	52	37	52	52	37	37	36
	ɛ_S=10^−5	95.62 s	33	50	34	50	50	34	34	33
Case Study 2	ɛ_S=10^−3	92.48 s	28	37	29	37	37	29	29	28
	ɛ_S=10^−4	95.90 s	33	50	34	50	50	34	34	33
	ɛ_S=10^−6	98.41 s	34	48	35	48	48	35	35	34
Case Study 3	A=0.01	106.51 s	36	52	37	52	52	37	37	36
	A=0.10	102.08 s	34	51	35	51	51	35	35	34
	A=1.00	108.50 s	35	53	36	53	53	36	36	35

Conclusions

In this paper, the formulation and numerical resolution of robust optimal control problems for compartmental epidemic models using a spectral approach have been presented. The proposed methodology has been described assuming the presence of mixed state-control constraints in the optimal control problem. However, it could be extended to more general formulations with arbitrary algebraic constraints.

To illustrate the effectiveness of the proposed approach, a mathematical model of COVID-19 transmission that includes a vaccination strategy has been considered, in which uncertainties involving the replication factor and the isolation rates of symptomatic and asymptomatic individuals have been assumed in an effort to find a robust optimal vaccination strategy.

Three solutions to the robust optimal control problem have been computed using three different values of the upper bound on the standard deviation of the proportion of susceptible individuals. In these numerical experiments several aspects of the solutions have been studied. First, the percentages of symptomatic individuals in a single day has been estimated to assess if they place undue stress on the healthcare system. Then, the proportion of susceptible individuals that would become immune through vaccination over a certain time period has been estimated to determine if herd immunity could be reached or not in this period. The probability that, at a given time instant, the percentage of symptomatic individuals exceeds a given threshold has been also calculated. In addition, a multivariate statistical analysis has also been conducted, more specifically, the correlation matrix for the distributions of symptomatic individuals at given time instants has been calculated. Furthermore, conditional probabilities have been estimated, in particular, the probability that the percentage of symptomatic individuals at a given time instant exceeds a given threshold assuming that this percentage exceeded another threshold at a preceding time instant. Finally, the Sobol’ indices have been calculated to understand which random parameters have a greater influence on the variability of the state variables at a given time instant.

Additional numerical experiments have been conducted in order to assess the effects of a decrease in the rate of identification and isolation of asymptomatic individuals on the behavior of the optimal control and state variables. Moreover, the sensitivity of the robust optimal control problem to the weighting parameter of the objective functional has also been analyzed.

The possibility of imposing an upper bound on the standard deviation of the state variables to reduce their variability, and the availability of several analysis tools, demonstrate that this methodology is a valuable tool to control and mitigate the effects of infectious diseases.

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.