Introduction

Mathematical models of infectious diseases have been extensively studied over the past decades [1,2,3,4,5,6,7,8,9,10,11]. The outbreak of infectious disease has caused mortality of millions of people as well as represents a noteworthy risk to public health and may cause heavy economic and social losses [12]. Needless to say, the quick spread of the coronavirus in 2020 throughout the world is a clear evidence of the importance of research on this topic [13, 14]. Among the several models of infectious diseases, the SIR (susceptible, infected, recovered) model, introduced in [15], is probably the most studied. In spite of the fact that, it appears to be basic, the SIR model has been explored in numerous viewpoints, for example, nonlinear examination of wave going of SIR pandemic model [16], effects of discretization schemes [17, 18], HIV on small scale population [19] variable population [20, 21], and optimal control [22,23,24,25,26,27] to cite a few examples.

The SIR model has been important to clarify several key issues. One of them is the fact that there is a threshold for the strength of disease spread, wherein there is an endemic behaviour of an infectious disease. Besides that this threshold can be changed by means of vaccination [28]. More specifically, by vaccinating a sufficient percentage of the population, the disease can be eradicated. In fact, many works have addressed with some success this issue using a constant vaccination rate [3, 29, 30]. In fact, in several public health practices, a constant value of vaccination rate is adopted as target during several public campaigns [31]. Recently, the author of [32] has applied optimization techniques to find suitable vaccination strategies in an age structured SIR model.

Due to economic pressure over the government agencies all around the world, suitable vaccination campaigns are difficult to implement [33]. This fact has made several researchers turn their attention to the optimal control problem [34]. The author of [22] has considered various deterministic optimal control models for SIR-epidemics. The author has used vaccination, quarantine, screening or health promotion campaigns as forms of control. In [35] the authors have applied an optimal control strategy during flu season, the SIR model was also used in that work.

In [36], the goal was to figure ideal vaccination designs for a quickly spreading disease in an urbanized profoundly mobile population, using different types of structured SIR models. In [37] the authors have developed a time-delay SIR model to consider that infected individuals can remain sick for a period of time defined over an interval with a lower bound. An optimal impulse control strategy has been employed, the cost function considered was composed by a cost associated with the state variables and a cost associated with the vaccination (that depends on the control effort and also on delayed states). In [12] the stability and the optimal control problem for the SIR model has been investigated. The objective functional to be minimized depends on the number of susceptible and infected individuals in addition to the control effort. The authors in [38] have considered the optimal control problem for the standard SIR model, the main target was the minimization of the infection eradication time. The authors of [26] present an optimal control approach for the SIR model based on a cost index that depends on the state variables, on the control effort that is weighted by a nonlinear function and it also depends on mixed terms (product between number of infected individuals and control effort). Optimal control has also been considered to deal with the SIR model under the assumption of imprecise parameters described as intervals [39]. In that case a linear cost function of infected individuals and control effort has been employed. One may cite also the work proposed in [40] that deals with the optimal control problem for spatiotemporal SIR models. The system was described by a set of partial differential equations and the considered control variable was the spatial and temporal distribution of the vaccine.

In the aforementioned articles there are important results for the research on epidemiology. Nevertheless, the interpretation of an optimal control law for an epidemic system has been little investigated. This seems an issue that has not received attention in the literature. In this paper, we have applied optimal control to eradicate infectious diseases described by SIR models. We also exploit the optimal control law to provide qualitative insights for control policies of vaccination. Using Pontryagin’s maximum principle [34] an optimal control law that minimizes the number of infected individuals and the vaccination rate has been obtained. However, there is not always a prior knowledge of the number of susceptible, infected, and recovered individuals required to formulate and solve the optimal control problem. In a model-free scenario, a strategy based on the analytic function is proposed, where prior knowledge of the scenario is not necessary. We have noticed that a very high level of vaccination rate in the beginning of control policy, with a systematic reduction in its level can be a general principle to be used. This principle is derived from the optimal control law. The considered approach is capable of eradicating the disease faster than a constant vaccination control method. The control employed in this paper can give some insights about how to act to eradicate infectious diseases in real life. Nowadays, the unique method to eradicate or mitigate the COVID-19 is the physical distancing measurements [11, 41]. After a development of a vaccine, there is no guarantee of worldwide cover due to economical and cultural restrictions [42]. Thus, it is still important to address such a problem for the COVID-19 pandemic. The results found here can be useful to stress the importance of a global and a general vaccination campaign. It also stresses the relevance of national and international cooperation to improve our preparedness for fighting pandemics such as COVID-19 [43].

The remainder of the paper is organized as follows. “Vaccination in the SIR Model ”shows the impact of vaccination on the dynamical behaviour of SIR. “Optimal Control” focuses on the optimal control problem. After that, in “Numerical Results” we solve numerically the system optimization problem. Thereafter, we display some results and compare with the solutions and cost obtained for a constant vaccination campaign and for the case without vaccination. Finally, in the last Section we conclude by discussing the main results of this paper.

Vaccination in the SIR Model

For a constant population, a normalized version of SIR model with vaccination can be seen as

$$\begin{aligned} \frac{ds}{dt}& = \mu (1-v) - \mu s - \beta is \nonumber \\ \frac{di}{dt}& = \beta is - \gamma i -\mu i, \end{aligned}$$
(1)

where \(\beta\) is the contact rate – also called transmission coefficient, it is the average number of contacts of an individual per unit of time, causing the transmission of the disease; \(\gamma\) is the recovery rate – it represents the rate of infected individuals per unit of time that passes to the recovered class; \(\mu\) is the renewal rate – it is the number of individuals who dies per unit of time and, in the same number, are born other susceptible individuals; v is vaccination rate; s is the proportion of susceptible individuals; i is the proportion of infected people.

When \(v =0\), Eq. (1) is equivalent to the SIR model [1]. The stability analysis of Eq. (1) has been conducted in [42]. Let be \(\beta _v=\beta (1-v)\). In order to guarantee the eradication of disease, we must hold

$$\begin{aligned} v > v_\mathrm {c} = 1- \frac{\gamma + \mu }{(\beta )}. \end{aligned}$$
(2)

Therefore, when \(v>v_\mathrm {c}\) the SIR model (1) presents eradication of infected individuals in steady-state regime. Other approaches create a new class, usually denoted as V, standing for vaccinated individuals [44]. For more details on the stability analysis of this model, please refer to [42].

Optimal Control

The optimal control problem consists of minimizing the objective functional:

$$\begin{aligned} J(I,v,t)= \int _{t_0}^{t_f} f_\mathrm {c}(I(\tau ),v(\tau ))d\tau , \end{aligned}$$
(3)

subject to the SIR model with vaccination

$$\begin{aligned} \frac{dS}{dt}& = \mu N(1-v) - \mu S - \beta I S/N, \nonumber \\ \frac{dI}{dt}& = \beta I S/N - \gamma I -\mu I, \nonumber \\ \frac{dR}{dt}& = \gamma I -\mu R+ v \mu N, \end{aligned}$$
(4)

where \(f_\mathrm {c}(\cdot )\) is a cost function; \(t_0\) is the initial time and \(t_f\) is the final time. In this case, we have considered the cost of the number of infected individuals I and the vaccination v. The description of parameters in Eq. (4) can be seen in [42].

The number of individuals with the disease (I) and the vaccination rate (v) to compose the cost function were chosen because they represent, respectively, the most important variable and in the majority of cases, most of the involved costs. Firstly, a minimal amount of investment required for vaccination is desired and also a minimal number of infected individuals. Second this allows one to calculate the cost function, with appropriated weights, \(\phi _1\) and \(\phi _2\). These weights can incorporate specific information concerning cost of vaccinations, vaccine campaigns and cost associated to each infected individuals. Thus, the cost function may be expressed as

$$\begin{aligned} f_c \left( I(\tau ), v(\tau ) \right) = \left( \phi _1 v(\tau )^2+ \phi _2I(\tau )^2 \right) . \end{aligned}$$
(5)

The basic quadratic form was chosen because, apart from its simplicity, it yields important features to the problem. Pontryagin’s maximum principle [45], a variational calculus method, has been employed in order to solve the optimal control problem, and we adopted the final time case. In the present case, the Hamiltonian is given by

$$\begin{aligned} \mathcal {H}_{(S,I,R,v,t)}& = \phi _1 v^2 + \phi _2 I^{2}+\left[ p _{1}\quad p _{2}\quad p _{3}\right] [S~I~R]^{\scriptscriptstyle \mathrm T } \nonumber \\& = \phi _1 v^2 + \nonumber \\&+ \phi _2 I^{2}+ \nonumber \\&+ p_1 [ \mu N(1-v) - \mu S - \beta I S/N ] + \nonumber \\&+ p_2 [ \beta I S/N - \gamma I -\mu I] \nonumber \\&+ p_3 [ \gamma I -\mu R+ v \mu N ], \end{aligned}$$
(6)

where \(p _{1}\), \(p _{2}\), \(p_{3}\) are the adjoint processes.

The control v has to minimize the Hamiltonian, so that, when the problem is unrestricted, it is required that

$$\begin{aligned} \left. \frac{\partial \mathcal {H}}{\partial v}\right| _{v=v^{*}}=0=2\phi _1 v -p_1 \mu N+p_3 \mu N. \end{aligned}$$
(7)

Rewriting (7) gives

$$\begin{aligned} v^*= \frac{\mu N(p_1 - p_3)}{2 \phi _1}. \end{aligned}$$
(8)

Note that Eq. (8) yields an analytic function for the optimal control law, where \(v^*\) depends directly on N and \(\mu\). The size of a population has already been investigated as an important parameter to indicate the persistence of infectious disease [46].

The actual determination of the optimal control, as shown in Eq. (8), requires the solution of an ordinary differential function, which is obtained by

$$\begin{aligned} -\frac{\partial \mathcal {H}}{\partial S}& = \frac{\mathrm {d} p_1}{\mathrm {d} t} = p_1 \mu + p_1 \beta I / N - p_2 \beta I /N \\ -\frac{\partial \mathcal {H}}{\partial I}& = \frac{\mathrm {d} p_2}{\mathrm {d} t}= \beta S / N(p_1 - p_2)+ p_2(\gamma +\mu ) - p_3 \gamma - 2 I \phi _1 \\ -\frac{\partial \mathcal {H}}{\partial R}& = \frac{\mathrm {d} p_3}{\mathrm {d} t} = p_3 \mu \end{aligned}$$

and with the boundary values

$$\begin{aligned} {[} S(t_0) \quad I(t_0) \quad R(t_0) ]& = [S_0 \quad I_0 \quad R_0] \end{aligned}$$
(9)
$$\begin{aligned} {[}p_1(t_f) \quad p_2(t_f) \quad p_3(t_f)]& = [ p_{\mathrm {1f}} \quad p_{\mathrm {2f}} \quad p_{\mathrm {3f}}], \end{aligned}$$
(10)

where \(S_0\), \(I_0\) and \(R_0\) are the state initial conditions. \(p_{\mathrm {1f}}\), \(p_{\mathrm {2f}}\) and \(p_{\mathrm {3f}}\) are the final conditions for the adjoint processes.

Numerical Results

In this section, we show numerical simulations of the optimal control applied to the SIR model. To compute the TPBVP (Two Point Boundary Value Problem) solution we have used the script bvp5c for Matlab [47]. Corresponding numerical performance in Matlab are also undertaken by [36].

The parameters related to the SIR model employed in this work are presented in Table 1. An arbitrary unit of time (u.t.) is adopted. It means, the user can set up the parameter according to the disease or scale-time of interest. For typical human diseases, this unit of time is usually adopted in days. We consider a simulation with 2 million individuals, with contact rate \(\beta = 0.08\), a recovery rate of \(\gamma = 1/24\). The critical value for vaccination \(v_\mathrm {c}\), using Eq. (2), is given by

$$\begin{aligned} v_\mathrm {c}=1-\frac{1/24+1/70}{0.08} = 0.3006. \end{aligned}$$

In the simulation of constant control, \(v_\mathrm {c}=0.40\) was used to guarantee the eradication of the disease. As an example, the worldwide measles vaccination coverage was 72% in 1998, which is considered insufficient to eradicate measles.

The choice of initial conditions represents a scenario where a small proportion, compared to the number of susceptibles, of infected individuals are present in a city with a population of two million. The adjoint processes were adjusted to minimize the maximum residual of the bvp5c algorithm.

Table 1 Parameters used in simulation
Full size table

Figures 1,2, 3 show the comparison of results for the number of infected individuals, vaccination and cost, respectively. In each figure, there are three situations: vaccination designed by optimal control, constant vaccination and without vaccination. The total cost in each situation was put in a per unit system. Figure 4 shows the adjoint process \(p_1\).

Fig. 1
figure 1

Number of infected individuals. The parameters of model and optimal control are presented in the Table 1. a Vaccination designed by optimal control (8). b Constant vaccination \(v_c = 0.40\). c Without vaccination. Notice that the same initial condition is adopted for all three cases. It is clear that the simulation performed with the optimal control presents the lowest peak of infected and the shortest time to eradicate the disease

Full size image
Fig. 2
figure 2

Comparison between vaccination in three situations. Optimal Control (black dashed line). Constant vaccination (dashed dotted red line). Without vaccination (dotted blue line). The constant vaccination rate is obtained by Eq. (2) and optimal control is obtained by Eq. (8). The y-axis is the vaccination rate and x-axis is time in an arbitrary unit. From this graph it is possible to notice an interesting behaviour of the optimal control. It starts with a very high vaccination rate, which allows a sharp reduction afterwards

Full size image
Fig. 3
figure 3

Comparison between cost obtained by (3) in three situations. Optimal control (black dashed line). Constant vaccination (dashed dotted red line). Without vaccination (dotted blue line). The total cost in each situation was put in a per unit system. The y-axis is the cost and x-axis is time in an arbitrary unit. In the beginning the cost of optimal control is higher, which is clearly related to an adoption of a much higher vaccination rate. However, after some time, this cost becomes lower than the cost provided by constant vaccination

Full size image
Fig. 4
figure 4

Adjoint process \(p_1\). This behaviour is close related to the optimal control seen in Fig. 2

Full size image

Another important feature of the application of the optimal control can be concluded from Fig. 5. Using the parameters of Table 1 and the stability analysis for the fixed point \(P_1\), an endemic state should occur when \(\beta >\gamma +\mu =0.056\). In this state, \(P_1\) is not stable and \(I>0\). We chose \(\beta\) varying from 0.056 to 0.1. The initial condition of the optimal control \(v^*(0)\) and the corresponding \(v_\mathrm {c}\) for each value of \(\beta\) are shown in Fig. 5. The initial value of the optimal control is higher than the critical value of vaccination when \(\beta\) increases. As the value of \(\beta\) is slightly greater than the minimum for many infectious diseases (see [3]), this fact can reveal an important practical orientation. The cost of vaccination campaigns can be optimized if the initial value is taken higher than the critical value for vaccination \(v_\mathrm {c}\).

Fig. 5
figure 5

Initial value of vaccination rate using the optimal control (–) and constant vaccination (-o-) for different values of \(\beta\). Notice that for \(\beta >0.06\), the vaccination rate designed by optimal control is higher than the critical value Eq. (2). It is also important to stress that in the majority of the cases the initial value is much higher. Moreover, values above 1 should be obviously saturated in 1 (maximum population cover of 100%)

Full size image

Vaccination Policy Based on Optimal Control Law

As presented, optimal control is able to provide satisfactory results. However, there is not always a prior knowledge of the number of susceptible, infected, and recovered individuals required to formulate and solve the optimal control problem. Therefore, constant vaccination has been used as a way to mitigate the effects of lack of knowledge of system states, but this approach produces a very high cost when compared to the cost of applying the optimal control.

To overcome the highlighted difficulties, a vaccination policy is proposed in which prior knowledge of the states is not necessary. The proposed approach starts with a high percentage of vaccination, to emulate the optimal control action. This percentage is then reduced to half of the initial value. From this moment on pulses of vaccination with constant interval and decreasing intensity are applied until the percentage of vaccination is zero. It is important to state that if the number of infected individuals does not reach zero, a small rate of vaccination can be continued applied.

Figures 6 and 7 show the comparison of results for vaccination and cost, respectively. In each figure, there are the situations: vaccination designed by optimal control, constant vaccination, without vaccination and proposed vaccination. The total cost in each situation was put in a per-unit system.

The weight for the number of infected individuals was varied and the results are summarized in Table 2. Optimal vaccination represents the minimum cost, except for \(\phi _2 = 0.1\), when the system without vaccination represents the minimum value. As can be seen, the final cost obtained by the proposed methodology is higher than the cost obtained by the optimal control and lower than the cost obtained by constant vaccination. It is noteworthy that the states were not used to formulate the control policy, which makes the scenario more realistic in view of the lack of accurate information about the system.

An important aspect that deserves our attention is related to the possibility of a second infection. It has been reported that it is likely that people can be reinfected with COVID-19. However, the degree of protective immunity conferred by infection by COVID-19 is currently unknown [48]. It means that at this stage of scientific development it would be very hard to incorporate such variables in the dynamics of the system. Nevertheless, one of the great advantages of compartmental models is the ability to easily incorporate new state variables. As soon as it would be possible to recognise the dynamics of the system, the loss of immunity can be included. This can be done, for instance, by adapting the SIR epidemic model with partial temporary immunity, as presented in [49].

Fig. 6
figure 6

Comparison between vaccinations. Proposed control policy (full green line). Optimal Control (black dashed line). Constant vaccination (dashed dotted red line). Without vaccination (dotted blue line). The constant vaccination rate is obtained by Eq. (2) and optimal control is obtained by Eq. (8). The proposed control policy (full green line) is particularly interesting for situations where the model is unknown or there is a substantial amount of uncertainty in the data to produce a blax-box model

Full size image
Fig. 7
figure 7

Comparison between cost obtained by (3). Proposed control policy (full green line). Optimal control (black dashed line). Constant vaccination (dashed dotted red line). Without vaccination (dotted blue line). The total cost in each situation was put in a per unit system. The y-axis is the cost and x-axis is time in an arbitrary unit. The performance of the Optimal Control is the best. However, it requires a precise model. Thus, the proposed policy can be an useful alternative in a model-free scenario

Full size image
Table 2 Cost for different values of \(\phi _2\) in four different situations for vaccination: optimal, proposed in this work, constant, and without vaccination at all
Full size table

Conclusions

This paper has presented an optimal control law to minimize the number of infected individuals and the vaccination rate. Furthermore, the strategy considered is based on an analytic function. Some aspects have been investigated. First, the number of infected individuals is analysed. The optimal vaccination campaigns present significant differences (Fig. 1). Both methodologies guarantee eradication of the disease, but with different peaks and eradication time. Optimal control takes one third of time to eradicate the disease when compared to constant vaccination. The maximum number of infected in the case of optimal control is around 50% lower than in the case of constant vaccination. The system without vaccination stabilizes around \(1.5 \times 10^{5}\) infected individuals.

Secondly, the vaccination for the three situations are quite different, as shown in Fig. 2. The optimal vaccination starts at a higher level than the constant vaccination and decreases smoothly until it reaches zero at \(t=300\). We believe this can indicate an important principle to design vaccination campaigns. The campaign should start with a high value and then it might be reduced throughout the time (see Figs. 2 and 5). Some similar results have been reached in [50], where the authors show the importance of a great investment in the start of a campaign. Thus, a great investment in vaccination campaigns reduces the cumulative cost in the long term. From this and the fact that there is not always a prior knowledge of the number of susceptible, infected and recovered individuals required to formulate and solve the problem of optimal control, a vaccination policy was proposed in which the prior knowledge of the states is not necessary. The proposed approach starts with a high percentage of vaccination, to emulate the optimal control action. This percentage is then reduced to half of the initial value. From this moment on pulses of vaccination with constant interval and decreasing intensity are applied until the percentage of vaccination is zero (see Fig. 6). Figure 7 presents the cumulative cost for the four studied cases. The results indicate that the absence of vaccination has the highest cost in the long term. In the short-term, the total cost of the optimal control is higher than the other cases. Besides that, the final cost obtained by the proposed methodology is higher than the cost obtained by the optimal control and lower than the cost obtained by constant vaccination, noting that neither prior knowledge of the number of susceptible, infected and recovered individuals was used in this approach.

Finally, another point is that the use of Pontryagin’s maximum principle yields an analytic function for the optimal control law. Thus, this conclusion can be expanded for a more realistic scenario the persistence of an infectious disease can also be described because of this relation. Although, in general there is no optimal design for vaccination campaigns, it is clear that economic pressures can be seen as a cost function, which make governmental agencies optimize their resources, and therefore, the vaccination campaign. Another interesting possibility to investigate is related to incorporating uncertainty from data and from computer simulation using interval arithmetic approaches for nonlinear dynamical systems as in [51,52,53].

Finally, it is important to stress that under the enormous pressure that global society has been facing to fight the COVID-19 pandemic, the results pointed here are based on the optimistic perspective for the development of a vaccine [54]. There is one lesson that, unfortunately, almost all countries have not learned: be prepared for a pandemic, as has been warned by many scientists and global leaders [55]. Therefore, it is really important to start to investigate how to deal with a massive and global vaccine campaign. The results here pointed out the importance of a robust response of vaccination to a quick eradication of the COVID-19. Future works should address other compartmental models, such as SEIR [41], and specific parameters for each country or region of the globe. It would be of great interest to incorporate other techniques to forecast COVID-19 as seen in [56].