Advanced modeling and parameter estimation algorithms form a solid background for the design of optimal strategies to control infectious diseases, which reduces illness and mortality rates. Vaccination, isolation, and public health education are examples of important control techniques that protect people at risk and make effective use of healthcare resources [1,2,3].

Timely control measures can mitigate the impact of outbreaks, prevent widespread transmission, and save lives. For instance, vaccination programs have been instrumental in controlling diseases such as measles, polio, and influenza [4,5]. Quarantine and isolation protocols were key in managing the spread of diseases like Ebola and COVID-19 [6]. Public health campaigns promoting handwashing and sanitation have significantly reduced the transmission of diseases such as cholera and dysentery. The eradication of smallpox is a prime example of how global vaccination campaigns can lead to complete elimination of a disease [7]. Similarly, the rapid response to the H1N1 influenza pandemic in 2009, including the development and distribution of vaccines, helped to control the spread of the virus and reduce its impact [8].

By analyzing data that use different models and algorithms, epidemiologists can forecast future incidence cases and evaluate various control strategies. Systematic preventive measures can help in reducing the spread of diseases. At first glance, choosing healthcare policies seems obvious, but in reality it is a very complicated task. One needs to put forward control strategies that are practical and, at the same time, have manageable consequences. At the onset of COVID-19, lockdowns were helpful, but they were not sustainable long term [9,10,11]. Thus, choosing the best intervention at the right time is critical [12,13].

The study in [14] introduced a two-stage epidemic model for the spread of COVID-19 and proposed optimal control strategies based on actual data and cost considerations. The research underscores the importance of contact tracing and isolation in minimizing the costs and effectively curbing the spread of a disease. Numerical simulations and model analysis provide actionable recommendations for public health authorities, highlighting the critical role of controlling the transmission rate in epidemic management.

The research in [15] modeled the spread of COVID-19 and assessed the impact of social intervention measures during the early outbreak phase, focusing on optimal control strategies and the identifiability of model parameters. It found that optimal control strategies, especially social distancing and self-isolation, significantly reduced transmission rates when implemented early. The study emphasized the importance of structural identifiability for accurate parameter estimation in COVID-19 models. It shows that implementing control measures effectively "flattens the curve" and lowers the burden on healthcare systems.

Another study, [16], focused on an $SIR$ model with saturated incidence and treatment rates, analyzing equilibrium points, bifurcation, and optimal control strategies that utilize vaccination and treatment with antiviral medication in order to contain the outbreak. The authors’ findings, derived from numerical simulations and efficiency analysis, demonstrated that vaccination control stands to reduce the cumulative number of infections more rapidly than control by antiviral treatment. This research underscored the value of mathematical modeling in epidemiology and the strategic implementation of vaccination to prevent disease transmission.

In the study of epidemic control, effective management of a disease spread is crucial, particularly at the onset of an outbreak. While the importance of vaccination in fighting infectious diseases is undeniable, it takes time to develop a vaccine for an emerging strain. Various parameters including environmental factors, immunity patterns, and behavioral changes impact the circulation of a virus. Social distancing and personal hygiene measures (non-medical interventions) play an important role in containing the disease at an early ascending stage. By optimizing the implementation of non-medical interventions over time the effectiveness of these interventions can be increased.

Another essential component of control and prevention is treatment with antiviral medications, which allows to reduce the period of infectiousness and/or reduce the disease fatality rate. To account for these two important types of control, we consider the following $SIR$ (Susceptible-Infectious-Removed) model [17] for early disease transmission: In system (1), we assume that individuals who have recovered from the disease gain immunity for the duration of the study and do not return to the susceptible class, $\mathcal{S}$. Additionally, we assume that natural birth and death rates balance one another and the removed class, $\mathcal{R}$, contains individuals who recovered from the disease and those who died. The two disease parameters, $\beta >0$ and $\gamma >0$ are transmission and recovery rates, respectively.

The focus of this research is on introducing optimal controls during initial weeks of a pandemic in order to delay and reduce the daily number of infections. This approach enables health centers and decision-making organizations to implement more effective operations. In what follows, we incorporate two different kinds of control in $SIR$ model [18] resulting in the system ${\displaystyle \frac{d\mathbf{x}}{dt}}=f(\mathbf{x},\mathbf{u})$, where Here $S\left(t\right):={\displaystyle \frac{\mathcal{S}\left(t\right)}{N}}$, $I\left(t\right):={\displaystyle \frac{\mathcal{I}\left(t\right)}{N}}$, and $R\left(t\right):={\displaystyle \frac{\mathcal{R}\left(t\right)}{N}}$ are normalized susceptible, infected and removed compartments, respectively, and N is the population of the region at the beginning of the study period. The function $u_1\left(t\right)$ represents nonmedical controls (social distancing, remote work, online education, restriction on travel, lockdowns, etc.), while $u_2\left(t\right)$ stands for treatment with antiviral medications and other medical interventions. A positive parameter, $\epsilon$, is the efficacy of antiviral treatments [19]. In the above, $\mathbf{x}:={[S,I,R]}^{\top }$, $\mathbf{u}:={[u_1,u_2]}^{\top }$, and the admissible set for each control function is In (2), the first control, $u_1\left(t\right)$, aims to change the disease transmission rate from $\beta$ to $\beta (1-u_1\left(t\right))$ while the second control, $u_2\left(t\right)$, is expected to reduce the period of infectiousness, which is ${\displaystyle \frac{1}{\gamma }}$ in the initial system (1). In combination, the two controls, $u_1\left(t\right)$ and $u_2\left(t\right)$, are meant to minimize the force of infection, $\beta (1-u_1\left(t\right))S\left(t\right)I\left(t\right)$, while keeping the costs at bay. The costs are considered in a general sense, which includes a negative impact on mental health, education, economy and on overall quality of life.

In Lemma 2.1 below we show that following the introduction of a time-dependent transmission rate, $\beta \left(t\right):=\beta (1-u_1\left(t\right))$, and a time-dependent recovery rate, $\gamma \left(t\right):=\gamma +\epsilon u_2\left(t\right)$, the model ${\displaystyle \frac{d\mathbf{x}}{dt}}=f(\mathbf{x},u)$ in (2) remains well-defined in the sense that the state variables, $S\left(t\right),I\left(t\right),R\left(t\right),$ originating in a positive octant do not leave the octant for all values of $t>0$. The proof of this lemma is similar to the argument in [20], where system (2) was considered with nonmedical controls only (that is, $u_2\left(t\right)=0$).

Lemma 2.1 [20]. Let $u_i\left(t\right),$$i=1,2$, be admissible control trajectories with $\mathbf{x}\left(t\right)$ satisfying ${\displaystyle \frac{d\mathbf{x}}{dt}}=f(\mathbf{x},\mathbf{u})$ given by (2) and the probability simplex in ${\mathbb{R}}^3$. Then $(S\left(t\right),I\left(t\right),R\left(t\right))\in {\Delta }^2$ for all $t\ge 0$.

Note that the argument in [20] implies that the conclusion of Lemma 2.1 is not limited to $\beta \left(t\right):=\beta (1-u_1\left(t\right))$ and $\gamma \left(t\right):=\gamma +\epsilon u_2\left(t\right)$. It is valid for any integrable nonnegative functions $\beta \left(t\right)$ and $\gamma \left(t\right)$. To optimize the implementation of both controls, $u_1\left(t\right)$ and $u_2\left(t\right)$, we consider the following objective functional According to system (2), one can integrate the first term to obtain where $C\left(\mathbf{u}\right):={[C_1\left(u_1\right),C_2\left(u_2\right)]}^{\top }$ is the assumed cost of control and $\lambda :={[{\lambda }_1,{\lambda }_2]}^{\top }$, ${\lambda }_1,{\lambda }_2>0$, is the regularization parameter (weight). As our numerical experiments show, the choice of the cost function, $C\left(\mathbf{u}\right)$, significantly influences the resulting control strategy. From practical standpoint, neither $u_1\left(t\right)$ nor $u_2\left(t\right)$ should take negative values. At the same time, the cost, $C_i\left(u_i\right)$, must increase dramatically as $u_i\left(t\right)$ approaches 1, the upper bound of the feasible set (3), since it is impossible to entirely eliminate the disease transmission $(u_1\left(t\right)=1)$. It is equally impossible to reach full capacity of antiviral treatment ($u_2\left(t\right)=1$) due to limitations of testing and other factors. Therefore, in our approach we impose the following assumptions on the cost functions, $C_1\left(u_1\right)$ and $C_2\left(u_2\right)$ [20]: These assumptions on the cost of control were first proposed in [20] for a special case when $u_2\left(t\right)=0$. The authors of [20] employed the techniques of machine learning to show that under assumptions (5), the global minimum of the Hamiltonian gives rise to the optimal control strategy, $u_1\left(t\right)$, which stays within the feasible set (3) for all values of $t\in [0,T]$. Assumptions (5) are the alternative to a more traditional cost function, $C\left(u\right)=u^2$, that is often used in control literature. However, $C\left(u\right)=u^2$ does not, generally, prevent the global minimum from becoming greater than 1 for some values of t, even for large weights, $\lambda .$

In this section we state and prove our main theoretical result.

Theorem 3.1. Let $u\in \mathcal{U}$ be an optimal control strategy with respect to the objective functional $J(\mathbf{x},\mathbf{u})$ defined in (4) and biological model $\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u})$, $\mathbf{x}\left(0\right)={\mathbf{x}}_{\mathbf{0}}$, introduced in (2), with $C\left(\mathbf{u}\right):={[C_1\left(u_1\right),C_2\left(u_2\right)]}^{\top }$ satisfying (5) and $\lambda :={[{\lambda }_1,{\lambda }_2]}^{\top }$, ${\lambda }_1,{\lambda }_2>0$. Then there is $\tau \in [0,T)$ such that for any $t\in (\tau ,T)$, the derivative, ${\displaystyle \frac{d{u}_i}{dt}}$, $i=1,2,$ exists and ${\displaystyle \frac{d{u}_i}{dt}}<0$. In other words, there is $\tau \in [0,T)$ such that for any $t\in (\tau ,T)$, both optimal controls, $u_1\left(t\right)$ and $u_2\left(t\right)$, are decreasing.

Proof. According to the Pontryagin’s Minimum Principle [21,22], if $u\in \mathcal{U}$ is an optimal control with respect to the objective functional $J(\mathbf{x},\mathbf{u})=h\left(\mathbf{x}\left(T\right)\right)+{\int }_0^{T}L(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right))dt$ and the system $\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u})$, $\mathbf{x}\left(0\right)={\mathbf{x}}_{\mathbf{0}}$, then there is a trajectory $\mathbf{p}\left(t\right)$ such that By the properties (5) of the cost, $C\left(\mathbf{u}\right)$, one has $C_i^{\prime }\left(u\right)>0$ for $u>0$ and ${lim}_{u\to 1^{-}}{C}_i\left(u\right)=\infty$, which prevents any $\mathbf{u}={[u_1,u_2]}^{\top }$, $u_i\ge 1$, $i=1,2,$ from being the optimal of $H(\mathbf{x},\mathbf{u},\mathbf{p})$ with respect to $\mathbf{u}$ at any time $t\in [0,T]$. Therefore, the Karush–Kuhn–Tucker (KKT) conditions for the optimal control problem (2), (3), (4) take the form As it follows from (2), (4), and (7), which yields To show that on some $(\tau ,T)$ the derivative, ${\displaystyle \frac{d{u}_1}{dt}}$, exists and ${\displaystyle \frac{d{u}_1}{dt}}<0$, we follow [20]. Conditions (K2) and (K3) imply that $(p_1-p_2)SI$ is differentiable and therefore continuous for any $t\in [0,T]$. From Lemma 2.1, one concludes that $S\left(t\right),I\left(t\right)>0$ as long as $S\left(0\right)$ and $I\left(0\right)$ are positive. On the other hand, since $p_1\left(T\right)=-1<0=p_2\left(T\right)$, there is ${\tau }_1\in [0,T)$ such that $p_1\left(t\right)-p_2\left(t\right)<0$ for all $t\in [{\tau }_1,T)$. Suppose at some point $t\in [{\tau }_1,T]$, the Lagrange multiplier, $q_1\left(t\right)$, is greater than zero. Then from (K4) it follows that $u_1\left(t\right)=0$. By (5), this implies ${\displaystyle \frac{d{c}_1}{d{u}_1}}\left(t\right)=0$, which means that in (14) ${\displaystyle \frac{d{c}_1}{d{u}_1}}\left(t\right)-q_1\left(t\right)<0$, while $-\beta (p_1-p_2)SI>0$. Hence, we arrive at the contradiction. Therefore, for any $t\in [{\tau }_1,T]$, one has $q_1\left(t\right)=0$ and ${\lambda }_1{\displaystyle \frac{d{c}_1}{d{u}_1}}=-\beta (p_1-p_2)SI$. By implicit function theorem, for $t\in ({\tau }_1,T)$ the derivative, ${\displaystyle \frac{d{u}_1}{dt}}$, exists and Taking into consideration (2), (4), and (7), one obtains Furthermore, from (4) one gets ${\partial }_{x}h\left(\mathbf{x}\right)={[-1,0,0]}^{\top }$. This together with (16) implies that $p_3\left(t\right)=0$ and the costate equations for $p_1\left(t\right)$ and $p_2\left(t\right)$ take the form: Combining (2) and (17), one can rewrite ${\left[S\left(t\right)I\left(t\right)(p_1\left(t\right)-p_2\left(t\right))\right]}^{\prime }$ as follows From (24) and (15), one concludes Since $p_1\left(T\right)=-1<0$ and $S\left(t\right),I\left(t\right)>0$, $\gamma +\epsilon u_2>0$, ${\lambda }_1{c}_1^{\prime \prime }\left(u_1\right)>0$ for $t\in [0,T]$, there exists ${\tau }_2\in [0,T)$ such that $p_1\left(t\right)<0$ and ${\displaystyle \frac{\beta p_1(\gamma +\epsilon u_2)SI}{{\lambda }_1{c}_1^{\prime \prime }\left(u_1\right)}}<0$ for all $t\in [{\tau }_2,T]$. Let $\tau =max({\tau }_1,{\tau }_2)$, then ${\displaystyle \frac{d{u}_1}{dt}}$ is negative in $(\tau ,T).$ As noted above, $p_3\left(t\right)=0$, therefore identity (14) yields Taking into account (17), one arrives at Integrating (21) from t to T and substituting $p\left(T\right)=0$, one gets As shown above, $p_1\left(t\right)-p_2\left(t\right)$ is negative on $[\tau ,T]$. Thus (3) and (22) imply that $p_2\left(t\right)>0$ for all $t\in [\tau ,T)$. Using the same argument as in the case of $u_1\left(t\right)$, one can now conclude that $q_2\left(t\right)$ in (20) is equal to zero on $[\tau ,T]$, that is, the constraint $u_2\left(t\right)\ge 0$ is not active for $t\in [\tau ,T)$ and ${\lambda }_2{\displaystyle \frac{d{c}_2}{d{u}_2}}=\epsilon p_{2}I$. By implicit function theorem, for $t\in (\tau ,T)$ the derivative, ${\displaystyle \frac{d{u}_2}{dt}}$, exists and From (2) and (17) one has since $p_1\left(t\right)<0$ for all $t\in [{\tau }_2,T]$ and $\tau \ge {\tau }_2$. This implies ${\displaystyle \frac{d{u}_2}{dt}}<0$ in $(\tau ,T)$, which completes the proof. □

Remark 3.2 By (5) and (7), ${\partial }_{\mathbf{u}}^{2}H(\mathbf{x},\mathbf{u},\mathbf{p})=\left[\begin{array}{cc}{\lambda }_1{c}_1^{\prime \prime }\left(u_1\right)& 0\\ 0& {\lambda }_2{c}_2^{\prime \prime }\left(u_2\right)\end{array}\right]$. Therefore, ${\partial }_{\mathbf{u}}^{2}H(\mathbf{x},\mathbf{u},\mathbf{p})$ is positive definite and $H(\mathbf{x},\mathbf{u},\mathbf{p})$ has a unique global minimum with respect to $\mathbf{u}$. From the proof of Theorem 3.1, it follows that both coordinates of the global minimum, $u_1\left(t\right)$ and $u_2\left(t\right)$, are guaranteed to be less than 1, pointwisely, but they are not guaranteed to be greater than 0 necessarily, which means the solution to our optimal control problem can be a local minimum rather than global. The reason the coordinates of the global minimum, $u_i\left(t\right)$, $i=1,2$, can potentially be less than zero for some values of t is that, counterintuitively, a smaller effective reproduction number, $r\left(t\right)$, in the SIR model does not always imply a smaller cumulative number of infected people, $S\left(0\right)-S\left(t\right)$. Hence, even though for system (2), the effective reproduction number, $r\left(t\right)=\beta (1-u_1\left(t\right))/(\gamma +\epsilon u_2\left(t\right))$, goes down with more control, it does not guarantee that $r\left(t\right)\ge \overline{r}\left(t\right)$ yields $S\left(t\right)\le \overline{S}\left(t\right)$ for every value of t. One can, however, show that if $r\left(t\right)\ge \overline{r}\left(t\right)$ and $r\left(t\right)$ is non-increasing, then $S\left(t\right)\le \overline{S}\left(t\right)$. This result is important in its own right. Its proof is given in Appendix.

Remark 3.3 Despite the fact that, theoretically, the solution to our optimal control problem can be a local minimum rather than global, in all numerical experiments presented in the next section, the optimal strategy is a unique global minimum. In other words, in all our experiments, the optimal control has been computed from the first-order necessary condition for unconstrained minimization, and nonnegativity constraints held without being enforced. For all cost functions, $C\left(\mathbf{u}\right(t\left)\right)$, satisfying (5), the global minimum of $H(\mathbf{x},\mathbf{u},\mathbf{p})$ with respect to $\mathbf{u}$ has nonnegative coordinates, $u_i\left(t\right)$, $i=1,2$. That is, $\mathbf{u}\left(t\right)={arg\; min}_{\mathbf{v}\in \mathcal{U}}H(\mathbf{x}\left(t\right),\mathbf{v}\left(t\right),\mathbf{p}\left(t\right))$ is equivalent to $\mathbf{u}\left(t\right)={arg\; min}_{\mathbf{v}\in {\mathcal{L}}^1[0,T]}H(\mathbf{x}\left(t\right),\mathbf{v}\left(t\right),\mathbf{p}\left(t\right))$ for all $t\in [0,T].$ This illustrates that conditions (5) lead to a practically justified mitigation scenario. Numerical simulations have also confirmed, as proven in Theorem 3.1, that both controls, $u_1\left(t\right)$ and $u_2\left(t\right)$, are decreasing toward the end of the study period.

In our numerical study, we use a second-order trust-region algorithm for nonlinear optimization ’lsqnonlin’ combined with ode15s built-in function to approximate the solution to optimal control problem (4) subject to compartmental model (2) and costate system (17). For every value of ${\mathbf{u}}_k$, we solve system (2) forward in time (starting with $\mathbf{x}\left(0\right)={\mathbf{x}}_{\mathbf{0}}$), to get ${\mathbf{x}}_k$ using ode15s. Then system (17) is solved back in time by ode15s to obtain ${\mathbf{p}}_k$. Given $({\mathbf{x}}_k,{\mathbf{u}}_k,{\mathbf{p}}_k)$, we find ${\mathbf{u}}_{k+1}$ by applying a second-order trust-region ’lsqnonlin’ algorithm.

Following [20], we consider three different cost functions, $C_{i,1}$, $C_{i,2}$, and $C_{i,3},$ satisfying conditions (5): In (26), the coefficients are chosen to minimize the weighted distance [20] The cost of control, $C_{i,1}\left(u\right)$, $C_{i,2}\left(u\right)$, and $C_{i,3}\left(u\right)$, is infinite as u approaches its ultimate value, 1. For comparison, we also use the cost function $C_{i,4}\left(u\right)=u^2$, for which (5) does not hold. The cost function $C_{i,4}\left(u\right)=u^2$ is popular in applications of control theory in epidemiology and other fields, since for this function the first order optimality condition is linear with respect to u. This is a useful property that simplifies numerical algorithms. However, the cost of control, $C_{i,4}\left(u\right)$, is finite at $u=1$, which is not realistic in real-world scenarios. Figure 4 and Figure 7 show that the global minimum, $\mathbf{u}\left(t\right)$, of the Hamiltonian, $H\left(\mathbf{x}\right(t),\mathbf{u}(t),\mathbf{p}(t\left)\right)$, does not stay in the range of $[0,1]$ if the cost is given by $C_{i,4}\left(u\right)=u^2$, especially for small values of ${\lambda }_i$, $i=1,2$. Thus, an explicit constraint $u_i\left(t\right)\le 1$ must be imposed in the case of $C_{i,4}\left(u\right)$. Even with the constraint $u_i\left(t\right)\le 1$, the optimal control function, $\mathbf{u}\left(t\right)$, often reaches the ultimate level [15], $u_i\left(t\right)=1$, which is not practical.

In all numerical experiments presented in this section, $C_{1,j}\left(u\right)=C_{2,j}\left(u\right)$, $j=1,2,3,4$. Therefore, moving forward, we omit the first index and set $C_{i,j}\left(u\right):=C_j\left(u\right)$. Comparison of the four cost functions in the interval $[-1,1]$ is illustrated in Figure 1.

In this study, numerical simulations are conducted for ${\lambda }_1$ and ${\lambda }_2$ equal to $0.1$, $0.05$, $0.01$, $0.001$ and ${10}^{-7}$. Three different scenarios are explored. First, there is only non-medical control, $u_1\left(t\right)$, (social distancing, behavioral changes, hand washing, etc.) in the system and treatment with antiviral medications, $u_2\left(t\right)$, is not available. Second, only control $u_2\left(t\right)$, treatment with antiviral medications, is applied; no social distancing. And third, controls $u_1\left(t\right)$ and $u_2\left(t\right)$, medical and non-medical, are used in combination. In our experiments, the population of the region, N, is assumed to be ${10}^7$. The initial number of infected individuals on day 1 is 200, and the duration of the study period is 120 days. Transmission rate, $\beta$, and recovery rate, $\gamma$, are assumed to be $0.3$, and $0.1$, respectively, leading to the basic reproduction number $r=3$. The efficacy of antiviral medication, $\epsilon$, is assumed to be $0.5$ when applicable.

To summarize, in this study, we investigated different control scenarios through theoretical analysis and numerical simulations. To account for two important types of control, social distancing and treatment with antiviral medications, the $SIR$ (Susceptible-Infectious-Removed) model [17] for an early ascending stage of an outbreak has been considered with the first control, $u_1\left(t\right)$, aimed to lower the disease transmission rate and the second control, $u_2\left(t\right)$, aimed to lower the period of infectiousness. In all experiments, the implementation of control strategies reduced the daily cumulative number of cases, $N-\mathcal{S}\left(t\right)$, and successfully "flattened the curve", $\mathcal{I}\left(t\right)$. The reduction in the cumulative cases was achieved by eliminating or delaying new cases. This delay is incredibly valuable as it provides public health organizations and researchers with more time to advance antiviral treatments and devise alternative preventive measures.

The main theoretical result of this paper, Theorem 3.1, concludes that the optimal control functions, $u_i\left(t\right)$, $i=1,2,$ may be increasing until some moment $\tau \in [0,T)$. However, for all $t\in [\tau ,T]$, the derivatives, ${\displaystyle \frac{d{u}_i}{dt}}$, become negative and $u_i\left(t\right)$ decline as t approaches T (possibly causing the number of newly infected people to grow). Numerical simulations presented in Section 4 confirm our theoretical findings. So, ideally, around the time $t=\tau$, the control strategy has to be upgraded and a vaccination campaign needs to start to ensure the epidemic wave does not rebound.