Optimal control analysis of COVID-19 stochastic model with comprehensive strategies

Abstract

In this paper, we construct a stochastic SVEAIR (Susceptible-Vaccinated-Exposed -Asymptomatic-Infected-Removed) epidemic model. It aimed to investigate COVID-19 control strategies through four different control methods, which consist of prophylaxis, vaccination control, rapid screening of people in exposure categories, and people identified as infected without screening. Firstly, we study the optimal control model and obtain the optimal control strategy using stochastic control theory. Secondly, the optimal control model with control term or without control term is numerically analyzed by using the implicit Runge-Kutta approximation method. At last, The numerical simulations verify the theoretical results and it show that control could help reduce infection and improve population health.

Keywords:

• Stochastic optimal control
• COVID-19
• stochastic SVEAIR model

1. Introduction

COVID-19 infection affected most countries and their economies in the world. Although many preventive mechanisms and other control measures have been put in place to reduce the spread of the disease, it was uncertain when this deadly infection would end in the population. Medical professionals, biologists and mathematicians were constantly working to identify effective vaccination, prevention and treatment measures to control coronavirus infection. Because there was so much variation in this infection, researchers were looking for some more effective vaccines to minimize infection.

Since its first discovery, COVID-19 has been a huge challenge, posing a serious threat to human life and health. Therefore, some strict measures must be taken to prevent virus transmission. On this basis, many scholars have applied mathematical models to analyze the epidemiological characteristics, transmission mechanisms and control strategies of COVID-19. The use of deterministic optimal control theory in epidemiological models dated back to the early days of the theory, see Marton and Wickwire (1974). The literature was extensive, ranging from simple compartment models, such as SIS models, to models involving distribution dynamics (heterogeneity in space, age, exposure rates, etc.). We mentioned (Behnke, 2000), where treatment and prevention (vaccination) was used as control in the SIR Model, Chen and Sun (2014a, 2014b), where control theory of treatment and vaccination was applied to the SIRS model in the network. There have been many research results on mathematical modelling and optimal control of the dynamics of COVID-19. For example, the dynamics of COVID-19 infection was studied through the mathematical model considered by the authors in Ullah and Khan (2020), in which the authors implement different combinations of controls to determine the best possible way to minimize infection. In Sun et al. (2020), the transmission trend of coronavirus infection in China and its modelling and prediction were discussed. In Khan et al. (2020), the authors described the dynamics of COVID-19 infection using isolation. In Asamoah et al. (2020), the mathematical model of severe acute respiratory syndrome coronavirus type 2 was studied through optimal control analysis. The application of an optimal control model using South African cases to coronavirus infection was studied (Gatyeni et al., 2022). There was still a lot of interesting work on coronavirus infection and its dynamic analysis and infection (Araz, 2021; Asamoah et al., 2021; Atangana, 2020; Razzaq et al., 2021).

The impact of environmental noise on health was increasingly attracting attention, for example, the transmission of dengue fever virus from mosquitoes to humans involved many random factors. The transmission mode of dengue fever infection was influenced by various random factors such as climate, migration, and changes in human behaviour. In order to consider these impacts in the study, researchers (Din et al., 2021; Khan et al., 2020; Qureshi & Atangana, 2019) developed many mathematical models to describe the dynamics of dengue fever through random disturbances involving human and mosquito populations. The study of the random COVID-19 model has similar aspects to the impact of random factors on dengue fever infection. In practice, environmental noise inevitably affected epidemics. More and more researchers were paying attention to this when studying mathematical models. Generally speaking, they used random models to study the impact of environmental noise. There were also some works on static optimization of stochastic epidemiological models in the literature. Several papers have analyzed vaccination in random environments, see Ball et al. (2004). In Kovacevic (2018), the authors analyzed the cost optimal long-term treatment and its stationary distribution of two stochastic SIS models. Recent work (Krause et al., 2018) numerically compared several fixed control strategies for networked SIS processes. However, the optimality of stochastic control problems in epidemiological models has so far seemed to be a largely neglected topic.

The main goal of this paper was to apply optimal control theory to the stochastic COVID-19 model. We considered four different control measures to minimise coronavirus infection and its further spread in the community. In the cost function, the nonlinear effect of control was considered. The optimal control path was characterized analytically by the stochastic Pontryagin's maximum principle. The numerical results indicated that under the proposed joint control, the disease transmission rate was low and the number of infections was minimized. The main contribution of the authors is to study COVID-19 using stochastic optimal control techniques. Its most important utility lies in maintaining a subtle balance between minimizing expenses and curbing the spread of the virus.

The rest of this paper is organized as follows. In the following section, we introduce a stochastic COVID-19 model with six compartments. We formulate the corresponding optimal control problem in Section 3. In Section 4, we carry out numerical experiments on the model system and discuss the experimental results.

2. Formulations of the model

Through the above detailed investigation of coronavirus infection research papers, we established a new stochastic dynamic mathematical model of coronavirus infection transmission. In order to establish the model of COVID-19, we divided the population by $N\left(t\right)$ into six different compartments, namely, susceptible $S\left(t\right)$, vaccinated $V\left(t\right)$, exposed $E\left(t\right)$, asymptomatic $A\left(t\right)$, infected $I\left(t\right)$ and recovered $R\left(t\right)$. The corresponding schematic diagram of the model dynamics was given in Figure 1. The deterministic model describing the dynamics of COVID-19 transmission in the community was represented by the following nonlinear differential equations (Shen, 2021): (1) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}}=\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu )S,\\ {\displaystyle \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}}=\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu )V,\\ {\displaystyle \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E,\\ {\displaystyle \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}}=\mathrm{\rho \delta E}-({\gamma }_1+\mu )A,\\ {\displaystyle \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}}=(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I,\\ {\displaystyle \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}}={\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}.\end{array}\end{array}$(1) where $\left(S\right(0),V(0),E(0),A(0),I(0),R(0)\in R_{+}^6$ was the initial value. In system (1(1) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}}=\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu )S,\\ {\displaystyle \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}}=\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu )V,\\ {\displaystyle \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E,\\ {\displaystyle \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}}=\mathrm{\rho \delta E}-({\gamma }_1+\mu )A,\\ {\displaystyle \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}}=(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I,\\ {\displaystyle \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}}={\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}.\end{array}\end{array}$(1) ), the growth rate of susceptible population was shown as Π. The parameter ${\text{β}}_i$ for i = 1, 2 described the rate of transmission of COVID-19 infection to susceptible individuals in asymptomatic and symptomatic patients, respectively, and assumed ${\text{β}}_1\ne {\text{β}}_2$. The exposed individuals who completed the incubation period were represented by the parameter δ, and became infected by joining the asymptomatic (with no clinical symptoms) class or symptomatic (with clinical symptoms) respectively, by $\mathrm{\rho \delta }$ and $(1-\rho )\delta$. The susceptible individuals were vaccinated at the rate of v, and the vaccine decay rate was ω. The parameter $0<\theta \le 1$ defined the vaccination effect of susceptible population on COVID-19 infection. Because it was difficult to identify the person infected with COVID-19, and the recovery of asymptomatic people was possible. We considered asymptomatic population recovery shown in the model by ${\gamma }_1$. Symptomatic infected patients died with infection rate d, and their recovery was expressed by ${\gamma }_2$. Individuals in each compartment in the model died naturally at the rate of μ.

Figure 1. Schematic diagram of the COVID-19 model, The meaning of the dashed arrow is as follows: $E\stackrel{\mathrm{\rho \delta E}}{⟶}A,E\stackrel{(1-\rho )\mathrm{\delta E}}{⟶}I,A\stackrel{{\text{β}}_1\frac{\mathrm{SA}}{N}+(1-\theta )\frac{{\text{β}}_1\mathrm{VA}}{N}}{⟶}E,I\stackrel{{\text{β}}_2\frac{\mathrm{SI}}{N}+(1-\theta )\frac{{\text{β}}_2\mathrm{VI}}{N}}{⟶}E$.

In the biological environment, it was more appropriate and reasonable to consider how environmental noise affects the spread of COVID-19. Numerous studies have shown that environmental fluctuations could affect the spread of epidemics in populations (see Bao et al., 2011; El Koufi et al., 2019; EL Koufi et al., 2022; Gao & Wang, 2019; Greenhalgh et al., 2016; Koufi et al., 2020; Zu et al., 2015). Thus, there were many ways to consider stochastic fluctuations in deterministic systems, one of which was to assume that an epidemic was affected by some small and standard stochastic fluctuations that could be represented by white noise. In this paper, we assumed that the COVID-19 model (1(1) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}}=\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu )S,\\ {\displaystyle \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}}=\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu )V,\\ {\displaystyle \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E,\\ {\displaystyle \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}}=\mathrm{\rho \delta E}-({\gamma }_1+\mu )A,\\ {\displaystyle \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}}=(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I,\\ {\displaystyle \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}}={\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}.\end{array}\end{array}$(1) ) was affected by environmental fluctuations, and the random disturbance was white noise type, which was proportional to the variable. We then introduced a stochastic SVEAIR model defined by the following stochastic differential equation system: (2) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu \left)S\right]\mathrm{d}t\\ +{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E]\mathrm{d}t+{\sigma }_3\mathrm{Ed}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta E}-({\gamma }_1+\mu )A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I]\mathrm{d}t+{\sigma }_5\mathrm{Id}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right),\end{array}\end{array}$(2) where $B_i\left(t\right)(i=1,2,3,4,5)$ was the standard Brownian motion defined on a complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathbb{P})$, filtered $\{{\mathcal{F}}_t{\}}_{t\ge 0}$ satisfying the usual conditions (i.e. it was increasing and right continuous, and ${\mathcal{F}}_0$ contained all $\mathbb{P}$-null sets), ${\sigma }_i(i=1,2,3,4,5,6)$ denoted the strength of the noise.

Theorem 2.1

For any initial value $\left(S\right(0),V(0),E(0),A(0),I(0),R(0)\in R_{+}^6$, there was a unique solution $(S^{\ast },V^{\ast },E^{\ast },A^{\ast },I^{\ast },R^{\ast })$ to model (2) on $t\ge 0$ and solution would remain in $R_{+}^6$ with probability one, namely $\left(S\right(t),V(t),E(t),A(t),I(t),R(t)\in R_{+}^6$ for all $t\ge 0$ almost surely.

By using the same steps of proof of Theorem 2.1 in EL Koufi et al. (2022) we could proved the above theorem. In the next section, we would focus on the optimal control problem of model (2(2) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu \left)S\right]\mathrm{d}t\\ +{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E]\mathrm{d}t+{\sigma }_3\mathrm{Ed}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta E}-({\gamma }_1+\mu )A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I]\mathrm{d}t+{\sigma }_5\mathrm{Id}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right),\end{array}\end{array}$(2) ), without discussing other properties.

3. Stochastic optimal control strategies

We applied the optimal control theory to the stochastic SVEAIR model given in system (2(2) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S-(v+\mu \left)S\right]\mathrm{d}t\\ +{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[\mathrm{vS}-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\delta +\mu )E]\mathrm{d}t+{\sigma }_3\mathrm{Ed}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta E}-({\gamma }_1+\mu )A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\mathrm{\delta E}-({\gamma }_2+\mu +d)I]\mathrm{d}t+{\sigma }_5\mathrm{Id}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right),\end{array}\end{array}$(2) ). We considered four different control measures. These controls could be defined as: Control variable $u_1$ was defined as prevention/isolation control, minimizing contact between healthy and infected people, regular hand washing, use of disinfectants and masks. Further, in areas with a large number of cases, avoid aggregation and limit their travel. The second control $u_2$ was vaccination control, where all possible individuals were vaccinated to further reduce transmission. Vaccination was the most effective way to reduce infection and further transmission, therefore it was recommended that individuals receive vaccines to reduce the risk of further transmission. The control $u_3$ denotes the rapid testing of the individuals in the exposed stage, and to identify the asymptomatic individuals. Upon identification through testing, the individuals should be isolated or restricted to their home etc to minimize the infection further. The control variable $u_4$ represents, the individuals that are not tested yet but identified the patients of COVID-19, either symptomatic or asymptomatic, can be treated and should be restricted to their places or hospitals, quarantines, etc. Therefore, $u_4$ should be considered as treatment control, while $b_1$ and $b_2$ should consider the rate at which asymptomatic and symptomatic individuals receive treatment. Through the above discussion, the following control system could be derived: (3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) The state variables were non-negative and the initial conditions were non-negative.

For the convenience of readers, we wrote (3(3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) ) as the following vector form, $\begin{array}{rl}x\left(t\right)& =\left[x_1\right(t),x_2(t),x_3(t),x_4(t),x_5(t),x_6(t){]}^{\prime },\\ u\left(t\right)& =\left[u_1\right(t),u_2(t),u_3(t),u_4(t){]}^{\prime }\end{array}$and $\begin{array}{r}\mathrm{d}x\left(t\right)=f\left(x\right(t),u(t\left)\right)\mathrm{d}t+g\left(x\right(t\left)\right)\mathrm{d}B\left(t\right),\end{array}$where the initial conditions was given as $x\left(0\right)=x_0$, f and g were vectors which comprises of the following components: $\begin{array}{rl}{f}_1\left(x\right(t),u(t\left)\right)& =\mathrm{\Pi }+\mathrm{\omega V}-\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}S(1-u_1)\\ & -(v{u}_2+\mu )S,\\ g_1\left(x\right(t),u(t\left)\right)& ={\sigma }_{1}S,\\ f_2\left(x\right(t),u(t\left)\right)& =v{u}_{2}S-(1-\theta )\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}V-(\omega +\mu )V,\\ g_2\left(x\right(t),u(t\left)\right)& ={\sigma }_{2}V,\\ f_3\left(x\right(t),u(t\left)\right)& =\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}S(1-u_1)\\ & +(1-\theta )\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}V\\ & -\left(\right(1-u_3)\delta +\mu )E,\\ g_3\left(x\right(t),u(t\left)\right)& ={\sigma }_{3}E,\\ f_4\left(x\right(t),u(t\left)\right)& =\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A,\\ g_4\left(x\right(t),u(t\left)\right)& ={\sigma }_{4}A,\\ f_5\left(x\right(t),u(t\left)\right)& =(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I,\\ g_5\left(x\right(t),u(t\left)\right)& ={\sigma }_{5}I,\\ f_6\left(x\right(t),u(t\left)\right)& ={\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R},\\ g_6\left(x\right(t),u(t\left)\right)& ={\sigma }_{6}R,\end{array}$The main purpose was to minimize the function, (4) $\begin{array}{rl}J(u_1,u_2,u_3,u_4)& =E{\int }_0^{T_f}\left(\frac{1}{2}\right[k_1{u}_1^2+k_2{u}_2^2+k_3{u}_3^2+k_4{u}_4^2]\\ & +k_{5}E+k_{6}A+k_{7}I\phantom{\frac{1}{2}})\mathrm{d}t,\end{array}$(4) subject to the control system (3(3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) ). The constant $k_i(i=1,\cdots ,7)$ described the weight or equilibrium constant, and $T_f$ was the final time. Because the control cost of intervention was nonlinear, the quadratic control function was realized. This means that there was no linear relationship between the cost of intervention and the effect of intervention. The controls defined above were Lebesgue integrable function and were bounded. We sought the optimal control $u_i^{\ast }(i=1,\cdot \cdot \cdot ,4)$ such that $\begin{array}{r}J\left(u_i^{\ast }\right)=\underset{U}{min}\left(J\right(u_i\left)\right),\end{array}$where U denoted to be the control set, (5) $\begin{array}{r}U=\{u_i:[0,T_f]\to [0,1],u_i\mathrm{is}\mathrm{Lebesgue}\mathrm{integrable}\}.\end{array}$(5) The existence of the stochastic optimal control could be obtained by using Yong and Zhou (1999). In proving the existence of optimal controls one typically sought a certain compactness structure. The demonstration of compactness depended on the specific control set used, usually piecewise continuous or bounded Lebesgue measurable functions. Subsequently, to show that the objective functional was bounded for all feasible controls, a maximized/minimized control sequence and a sequence of correlated states could be constructed and shown to converge in the appropriate feasible space. This principle transformed the control problem of finding the maximum objective functional satisfying the state SDE and the initial condition into the problem of optimizing the Hamiltonian function point by point. We defined the Hamiltonian function according to Yong and Zhou (1999). A necessary condition for $u_i$ could be constructed by maximizing the Hamiltonian at $u_i^{\ast }$, which was called the optimality condition. For the adjoint equation, we calculated the partial derivative of the state variable and took the final time condition of the adjoint variable, which was called the transversality conditions.

Theorem 3.1

The given optimal controls $u_i^{\ast }$ and solutions $S^{\ast },V^{\ast },E^{\ast },A^{\ast },I^{\ast },R^{\ast }$ of the control system (3(3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) ) that minimized $J\left(u_i^{\ast }\right)$ over U. Then there existed adjoint variables p statifying $\begin{array}{r}\mathrm{d}p=-\frac{\mathrm{\partial }H}{\mathrm{\partial }x}\mathrm{d}t+q\left(t\right)\mathrm{d}B\left(t\right)\end{array}$with the transversality conditions, $p\left(T_f\right)=0$.

The optimality condition was given by $\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_i}=0,i=1,\dots ,4$.

Furthermore, we had the controls (20(20) $\begin{array}{rl}& u_1^{\ast }=min\{max\{0,\frac{1}{k_1}(p_3^{\ast }-p_1^{\ast }\left){\displaystyle \frac{{\text{β}}_1{A}^{\ast }+{\text{β}}_2{I}^{\ast }}{N^{\ast }}}{S}^{\ast }\right\},1\},\end{array}$(20) )–(23(23) $\begin{array}{rl}& u_4^{\ast }=min\{max\{0,\frac{1}{k_4}(p_4^{\ast }{b}_1{A}^{\ast }+p_5^{\ast }{b}_2{I}^{\ast }\left)\right\},1\}.\end{array}$(23) ).

Proof.

Referring to the proof process of Yong and Zhou (1999), we briefly gave the following proof summary. The necessary condition of optimal control came from the stochastic Pontryagin's maximum principle. To use the stochastic maximum principle, we defined the Hamiltonian $H_m(x,u,p,q)$ as follows : (6) $\begin{array}{r}H(x,u,p,q)=〈f(x,u),p〉-L(x,u)+〈g\left(x\right),q〉,\end{array}$(6) where the symbol $〈\cdot ,\cdot 〉$ represented the Euclidean inner product, $p=[p_1,\cdot \cdot \cdot ,p_6]$ and $q=[q_1,\cdot \cdot \cdot ,q_6]$ were conjugate vectors. From the random maximum principle: (7) $\begin{array}{rl}& \mathrm{d}{x}^{\ast }\left(t\right)={\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }p}}\mathrm{d}t+g\left(x^{\ast }\right)\mathrm{d}B\left(t\right),\end{array}$(7) (8) $\begin{array}{rl}& \mathrm{d}{p}^{\ast }\left(t\right)=-{\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }x}}\mathrm{d}t+q\left(t\right)\mathrm{d}B\left(t\right),\end{array}$(8) (9) $\begin{array}{rl}& H(x^{\ast },u^{\ast },p,q)=\underset{u\in U}{max}H(x^{\ast },u,p,q),\end{array}$(9) where $x^{\ast }\left(t\right)$ was the optimal path of $x\left(t\right)$. The initial condition of Equation (7(7) $\begin{array}{rl}& \mathrm{d}{x}^{\ast }\left(t\right)={\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }p}}\mathrm{d}t+g\left(x^{\ast }\right)\mathrm{d}B\left(t\right),\end{array}$(7) ) was: (10) $\begin{array}{r}{x}^{\ast }\left(0\right)=x_0.\end{array}$(10) The terminal conditions of Equation (8(8) $\begin{array}{rl}& \mathrm{d}{p}^{\ast }\left(t\right)=-{\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }x}}\mathrm{d}t+q\left(t\right)\mathrm{d}B\left(t\right),\end{array}$(8) ) were: (11) $\begin{array}{r}p\left(T_f\right)=0.\end{array}$(11) Therefore, Equation (9(9) $\begin{array}{rl}& H(x^{\ast },u^{\ast },p,q)=\underset{u\in U}{max}H(x^{\ast },u,p,q),\end{array}$(9) ) showed that the optimal control $u^{\ast }\left(t\right)$ was a function of $p\left(t\right),q\left(t\right)$ and $x^{\ast }\left(t\right)$, namely, (12) $\begin{array}{r}{u}^{\ast }\left(t\right)=\varphi (x^{\ast },p,q)\end{array}$(12) where ϕ was determined by Equation (9(9) $\begin{array}{rl}& H(x^{\ast },u^{\ast },p,q)=\underset{u\in U}{max}H(x^{\ast },u,p,q),\end{array}$(9) ).

The maximum principle transformed the system (3(3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) ) and (4(4) $\begin{array}{rl}J(u_1,u_2,u_3,u_4)& =E{\int }_0^{T_f}\left(\frac{1}{2}\right[k_1{u}_1^2+k_2{u}_2^2+k_3{u}_3^2+k_4{u}_4^2]\\ & +k_{5}E+k_{6}A+k_{7}I\phantom{\frac{1}{2}})\mathrm{d}t,\end{array}$(4) ) into the pointwise minimization problem of the Hamiltonian H of the controller $u_i$. Therefore, the Hamiltonian function corresponding to Equations (8(8) $\begin{array}{rl}& \mathrm{d}{p}^{\ast }\left(t\right)=-{\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }x}}\mathrm{d}t+q\left(t\right)\mathrm{d}B\left(t\right),\end{array}$(8) ) and (9(9) $\begin{array}{rl}& H(x^{\ast },u^{\ast },p,q)=\underset{u\in U}{max}H(x^{\ast },u,p,q),\end{array}$(9) ) was defined as: $\begin{array}{rl}H& =\frac{1}{2}[k_1{u}_1^2+k_2{u}_2^2+k_3{u}_3^2+k_4{u}_4^2]+k_{5}E+k_{6}A+k_{7}I\\ & +p_1[\mathrm{\Pi }+\mathrm{\omega V}-\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}S(1-u_1)\\ & -(v{u}_2+\mu )S\phantom{\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}]+q_1{\sigma }_{1}S,\\ & +p_2[v{u}_{2}S-(1-\theta )\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}V-(\omega +\mu \left)V\right]\\ & +q_2{\sigma }_{2}V,\\ & +p_3\left[\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}S\right(1-u_1)\\ & +(1-\theta )\frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}V-\left(\right(1-u_3)\delta +\mu )E]\\ & +q_3{\sigma }_{3}E,\\ & +p_4[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]+q_4{\sigma }_{4}A,\\ & +p_5\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]+q_5{\sigma }_{5}I,\\ & +p_6[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]+q_6{\sigma }_{6}R.\end{array}$According to the random maximum principle, (13) $\begin{array}{r}\mathrm{d}{p}^{\ast }\left(t\right)=-{\displaystyle \frac{\mathrm{\partial }H(x^{\ast },u^{\ast },p,q)}{\mathrm{\partial }x}}\mathrm{d}t+q\left(t\right)\mathrm{d}B\left(t\right),\end{array}$(13) where (14) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_1}{\mathrm{d}t}}=(p_3-p_1)\frac{({\text{β}}_1{A}^{\ast }+{\text{β}}_2{I}^{\ast })}{N^{\ast }}(1-u_1)-p_1(v{u}_2+\mu )\\ & +q_1{\sigma }_1+p_{2}v{u}_2,\end{array}$(14) (15) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_2}{\mathrm{d}t}}=p_1\omega +(p_3-p_2)(1-\theta )\frac{({\text{β}}_1{A}^{\ast }+{\text{β}}_2{I}^{\ast })}{N^{\ast }}\\ & -p_2(\omega +\mu )+q_2{\sigma }_2,\end{array}$(15) (16) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_3}{\mathrm{d}t}}=k_5-p_3\left[\right((1-u_3)\delta +\mu \left)\right]+q_3{\sigma }_3\\ & +p_4\mathrm{\rho \delta }{u}_3+p_5(1-\rho )\delta u_3,\end{array}$(16) (17) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_4}{\mathrm{d}t}}=k_6-p_4({\gamma }_1+\mu +b_1{u}_4)+q_4{\sigma }_4+p_6{\gamma }_1,\end{array}$(17) (18) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_5}{\mathrm{d}t}}=k_7-p_5({\gamma }_2+\mu +d+b_2{u}_4)+q_5{\sigma }_5+p_6{\gamma }_2,\end{array}$(18) (19) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{p}_6}{\mathrm{d}t}}=-p_6\mu +q_6{\sigma }_6,\end{array}$(19) together with auxiliary initial and terminal conditions $\begin{array}{rl}& S^{\ast }\left(0\right)=S_0,V^{\ast }\left(0\right)=V_0,E^{\ast }\left(0\right)=E_0,A^{\ast }\left(0\right)=A_0,\\ & I^{\ast }\left(0\right)=I_0,R^{\ast }\left(0\right)=R_0,p_i\left(T_f\right)=0(i=1,\dots ,6).\end{array}$Now taking the differential of the Hamiltonian equation with respect to $u_i$ yielded the following optimal control $u_i^{\ast }(i=1,\dots ,4)$ (20) $\begin{array}{rl}& u_1^{\ast }=min\{max\{0,\frac{1}{k_1}(p_3^{\ast }-p_1^{\ast }\left){\displaystyle \frac{{\text{β}}_1{A}^{\ast }+{\text{β}}_2{I}^{\ast }}{N^{\ast }}}{S}^{\ast }\right\},1\},\end{array}$(20) (21) $\begin{array}{rl}& u_2^{\ast }=min\{max\{0,\frac{1}{k_2}(p_1^{\ast }-p_2^{\ast }\left)V^{\ast }{S}^{\ast }\right\},1\},\end{array}$(21) (22) $\begin{array}{rl}& u_3^{\ast }=min\{max\{0,\frac{1}{k_3}[(p_5^{\ast }-p_4^{\ast })\rho -p_3^{\ast }-p_5^{\ast }\left]\delta E^{\ast }\right\},1\},\end{array}$(22) (23) $\begin{array}{rl}& u_4^{\ast }=min\{max\{0,\frac{1}{k_4}(p_4^{\ast }{b}_1{A}^{\ast }+p_5^{\ast }{b}_2{I}^{\ast }\left)\right\},1\}.\end{array}$(23) This proof ended.

4. Numerical results

In this section, we performed numerical experimentations to analyze these results and also to see the effect of control interventions on disease dynamics. We numerically solved the control system (3(3) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}S\left(t\right)=[\mathrm{\Pi }+\mathrm{\omega V}-{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S(1-u_1)\\ -(v{u}_2+\mu )S\phantom{{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}}]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}V\left(t\right)=[v{u}_{2}S-(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-(\omega +\mu \left)V\right]\mathrm{d}t\\ +{\sigma }_{2}V\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}E\left(t\right)=\left[{\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}S\right(1-u_1)\\ +(1-\theta ){\displaystyle \frac{({\text{β}}_{1}A+{\text{β}}_{2}I)}{N}}V-\left(\right(1-u_3)\delta +\mu )E]\mathrm{d}t\\ +{\sigma }_{3}E\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}A\left(t\right)=[\mathrm{\rho \delta }{u}_{3}E-({\gamma }_1+\mu +b_1{u}_4)A]\mathrm{d}t+{\sigma }_{4}A\mathrm{d}{B}_4\left(t\right),\\ \mathrm{d}I\left(t\right)=\left[\right(1-\rho )\delta u_{3}E-({\gamma }_2+\mu +d+b_2{u}_4)I]\mathrm{d}t\\ +{\sigma }_{5}I\mathrm{d}{B}_5\left(t\right),\\ \mathrm{d}R\left(t\right)=[{\gamma }_{1}A+{\gamma }_{2}I-\mathrm{\mu R}]\mathrm{d}t+{\sigma }_{6}R\mathrm{d}{B}_6\left(t\right).\end{array}\end{array}$(3) )-(4(4) $\begin{array}{rl}J(u_1,u_2,u_3,u_4)& =E{\int }_0^{T_f}\left(\frac{1}{2}\right[k_1{u}_1^2+k_2{u}_2^2+k_3{u}_3^2+k_4{u}_4^2]\\ & +k_{5}E+k_{6}A+k_{7}I\phantom{\frac{1}{2}})\mathrm{d}t,\end{array}$(4) ). The values were as follows : $\mathrm{\Pi }=0.04$, $\mu =0.04$, ${\text{β}}_1=0.60$, ${\text{β}}_2=0.65$, $\delta =0.20$, ${\gamma }_1=0.50$, ${\gamma }_2=0.30$,d = 0.15, $\omega =0.02$, $\upsilon =0.08$, $\rho =0.50,\theta =0.1,b_1=0.3,b_2=0.4,k_1=0.25,k_2=1,k_3=0.3,k_4=1.5,k_5=0.6,k_6=0.1,k_7=\mathrm{0.2.}$ The initial value of the variable was taken as $N\left(0\right)=100$, $S\left(0\right)=57$, $V\left(0\right)=3$, $E\left(0\right)=28$, $I\left(0\right)=8$, $A\left(0\right)=2$, $R\left(0\right)=2$. All the numerical simulations are performed in MATLAB. The implicit Runge-Kutta approximation was used to obtain the graphical results shown in Figures 2–6. The left figure in Figure 2 showed the trajectory of the optimal prevention/isolation control curve, while the right figure showed the optimal trajectory curve of vaccine administration control. The left figure in Figure 3 showed the optimal detection control trajectory curve, while the right figure showed the optimal treatment control curve trajectory. Based on the optimal control curve trajectory in Figures 2-3, we simulated the dynamic trajectory curves of susceptible individuals, vaccinated individuals, exposed and asymptomatic individuals, infected individuals, and recovered individuals. The curve trajectory results showed that the control is very effective in eliminating diseases. In Figure 4, the left figure simulated the dynamic curve of the number of susceptible individuals, while the right figure showed the dynamic trajectory curve of the number of vaccinated individuals. In Figure 5, the left figure simulated the dynamic trajectory curve of the number of exposed individuals, while the right figure simulated the dynamic trajectory curve of the number of asymptomatic individuals. In Figure 6, the left figure simulated the dynamic trajectory curve of the number of infected individuals, while the right figure simulated the dynamic trajectory curve of the number of recovered individuals. From Figures 4–6, it could be seen that the control considered in this numerical simulation can effectively reduce the number of infected individuals and maximize the number of recovered individuals. Maintaining social distancing, wearing masks, regularly washing hands with disinfectants, vaccinating individuals, rapid testing, and possible early treatment could best reduce coronavirus infections in the community. From our random optimal control results, it could be seen that by implementing the above control, the number of infected cases decreases faster.

Figure 2. Optimal prevention/isolation control and vaccination control.

Figure 3. Optimal detection Control and treatment control.

Figure 4. Trajectory of changes in the number of susceptible individuals and vaccinated individuals.

Figure 5. Trajectory of changes in the number of exposed and asymptomatic individuals.

Figure 6. Trajectory of changes in the number of infected and recovered individuals.

5. Conclusion

We used stochastic optimal control theory to establish a mathematical model of coronavirus infection. Optimal control models had been implemented and formulated for controls that may be considered effective. The optimal control system was generated by four different control methods, and the desired results of the optimal control problem were obtained. The graphic results were obtained by combining the optimal control system, optimal control characteristics and adjoint equations. We propose some new results on the optimal control of random variables in the SVEAIR infectious disease model. The work is also dedicated to proposing various control strategies. These suggested techniques help us control the spread of diseases and save costs. By combining unpredictable factors such as timely detection and isolation of infectious sources, strengthening personal protection, and vaccination, the current work has been further extended. Some exciting and open topics are worth considering in the future. For example, using fractional order delayed differential equations to more finely characterize the process of disease propagation, and so on.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The work was supported by Yulin High-tech Zone Science and Technology Planning Project (CXY-2021-61, CXY-2021-65).