1 Introduction

On 31 December 2019, a novel coronavirus is formally reported in Wuhan, Hubei Province, China [1, 2]. On 31 January 2020, the World Health Organization formally knew this newly emerging virus as a global threat. In a very short time, this virus reached the most remote corner of the world. According to the WHO, human-to-human communication is the most common way of disease transmission. This fact highly affects the daily life of humans around the world.

This global pandemic leads to several economic losses and unemployment. For example, in [3, 4], the authors discussed effects of Covid-19 on the tourism industry. According to the global statistics, the tourism industry experienced the lowest demand in the last decade [5].The educational system was also affected severely by the Covid-19 pandemic [6,7,8]. Most schools and universities are experiencing virtual classrooms. In [9], the authors suggest that the virtual classroom trends remain even after the Covid-19 inhibition. The social and economic consequences of the Covid-19 pandemic have motivated researchers to study different methods to mitigate the spread of Covid-19.

The mathematical models of natural phenomena play a vital role for government and decision-makers. For example, mathematical models are employed to predict weather patterns [10, 11]. From initial days of the Covid-19 pandemic, scientists make a great effort to obtain the mathematical model of coronavirus [12, 13]. In [14,15,16], the author employs the classical SIR (susceptible–infectious–recovered) to describe the transmission dynamics of coronavirus in different societies. However, due to the complex behavior of coronaviruses, such as transmission power and asymptotic symptoms, this classical model is not able to reflect the behavior of the disease. In [17,18,19], the SEIR (susceptible–exposed–infectious–recovered) model is suggested to describe the transmission dynamics of Covid-19. In [20], the authors proposed the dynamical model of coronavirus in Mexico. The official data are employed to calibrate the model of disease, and the stability of the system is analyzed. In [21], the authors proposed a modified SEIR to model the disease outbreak in Italy and Spain for prediction and control purposes. The official data are employed to update the parameters of the modified SEIR model. In [22], the authors obtain the epidemiological model of the outbreaks in Brazil. The effect of social distance is analyzed in the mitigation of disease spread. In [23,24,25,26,27], the authors consider the effect of vaccination in modeling the disease outbreak. These models help us to understand how global vaccination reduces the rate of mortality.

Several researchers devoted their effort and attention to the parameter and state estimation of the epidemiological model. For example, in [28], the authors utilized neural networks and LSTM methods to predict the number of infected cases in India. In [29], the authors use the least square error method to estimate the unknown parameters of SIR/SQAIR (susceptible–quarantined–infected with no symptoms–infected with symptoms–recovered) models based on the actual data of South Korea. As a reliable method, the Kalman filter is widely employed to estimate the unknown parameters of the noisy processes. In [30, 31], the authors proposed a Kalman filter to estimate immeasurable parameters of epidemiological models. In [32], the authors try to predict the state variables of the epidemiological model of Covid-19 in Spain, France, and Italy based on the Auto-Regressive Integrated Moving Average method. In [33], based on the confirmed official data, the authors utilized deep learning methods to predict the spread of Covid-19 in six countries, including Australia, China, France, Italy, the USA, and Spain. The results of this research demonstrate the potential application of deep learning methods in the prediction of epidemiological model parameters. In [34], the fuzzy logic is proposed to predict the state variables of Covid-19 pandemics in the Indian society based on the official confirmed data.

In the last year, the design of a controller for the nonlinear epidemiological model of Covid-19 turns into a popular topic [35,36,37,38,39]. In [40], the authors proposed an optimal controller based on the Pontryagin principle to mitigate the spread of Covid-19. In [41], the authors develop a nonlinear model predictive control for the nonlinear epidemiological model of Covid-19 in Italy. The authors realized that social distance and restriction on mobility prevent the Covid-19 outbreak in the Italian society. In [42], the authors employ a Lyapunov-based method to investigate the effect of vaccination on the inhibition of Covid-19 in Canada. In [43], the authors design a controller for the SIR model based on the feedback linearization method. In this paper, the vaccination is considered as a control input signal.

According to the author’s best knowledge, the problem of robust sliding mode controller design based on the states and parameters estimation of the Covid-19 model is not studied yet. In this research, a novel mathematical model is introduced to reflect the behavior of Covid-19 in the Iranian and Russian communities. According to the actual data provided by the Iranian and Russian ministries of health, the transmission rates of the proposed model are estimated by the LSTM algorithm. In the LSTM method, the algorithm uses gates to control the memorizing process. Therefore, the problem of long-term dependency is prevented. In continuation, the immeasurable and noisy state variables of the proposed model are estimated by the EKF algorithm. To overcome the uncertainties and measurement noises in the nonlinear model of Covid-19, a robust sliding mode controller is suggested. Inspired by the above discussion, the main outcomes of this research are highlighted in below:

  • A novel mathematical model is proposed to describe the behavior of the Covid-19 pandemic.

  • Based on the actual data, the parameters of the model are estimated by the LSTM algorithm.

  • An extended Kalman Filter is suggested to estimate immeasurable state variables of model.

  • A robust sliding mode controller is proposed to mitigate the spread of Covid-19 through vaccination, social distance and facial masks, and medical treatment.

  • Global vaccination leads to herd immunity in the long term.

  • During the global vaccination, social distance and facial masks are obligatory.

The rest of this paper is organized as follows. In the second section, the mathematical model of Covid-19 with respect to the vaccination, social distance and facial masks, and medical treatment is introduced. In the third section, the parameter and state estimation are carried out by LSTM and EKF algorithms, respectively. In the next section, a robust sliding mode controller and Lyapunov stability analysis are developed. Finally, simulation results and concluding remarks are presented.

2 Mathematical model

In this section, the behavior of coronavirus under vaccination, social distance and facial masks, and medical treatment is expressed by mathematical language. In this epidemiological model, the population of Iranian and Russian societies are expressed with six states, including susceptible (S), infected (I), quarantined (Q), hospitalized (H), recovered (R), and finally extinct (E).

In Fig. 1, the transmission dynamics of the SIQHRE epidemiological model are presented. Normally, the RT-PCR test is the only method to detect coronavirus. The individuals with positive RT-PCR test move to the infected compartment with the rare of \(\alpha \).

Fig. 1
figure 1

The SIQHRE model

The infected cases with severe symptoms as pneumonia move to the hospitalized group at the rate of \(\gamma \). In this condition, the hospitalized cases that lose their life move to the extinct group with the rate of \(\varepsilon \); otherwise, they move to the recovered group with the rate of \(\phi \).

The infected cases with mild symptoms such as dry cough, tiredness, and fever move to the quarantined group with the rate of \(\lambda \). During this period, the quarantined cases with severe symptoms such as lung damage move to the hospitalized group at the rate of \(\theta \); otherwise, they move to the recovered group at the rate of \(\mu \).

In this research, the vaccination (u1), social distance and facial masks (u2), and medical treatment (u3) are considered as control input signals. Inspired by the above discussion, the mathematical model of the Covid-19 pandemic is presented by a set of six ordinary differential equations as follows

$$\begin{aligned} S(k+1)&=S(k)- \alpha I(k)S(k) \end{aligned}$$
(1)
$$\begin{aligned} I(k+1)&=I(k)+ \alpha I(k)S(k)-(\gamma +\lambda )I(k)+u2 \end{aligned}$$
(2)
$$\begin{aligned} Q(k+1)&=Q(k)+\lambda I(k)-(\theta +\mu )Q(k) \end{aligned}$$
(3)
$$\begin{aligned} H(k+1)&=H(k)+\gamma I(k)+\theta Q(k)-(\phi +\varepsilon )H(k) \end{aligned}$$
(4)
$$\begin{aligned} R(k+1)&=R(k)+\phi H(k)+\mu Q(k)+u1 \end{aligned}$$
(5)
$$\begin{aligned} E(k+1)&=E(k)+\varepsilon H(k)+u3 \end{aligned}$$
(6)

Remark 1

Since coronavirus statistics are updated on a daily basis, the SIQHRE model is expressed in the discrete-time framework with sampling period of one day.

Remark 2

It should be mentioned that all the state variables (populations) take positive values. According to the epidemiological dynamics, the population of community during Covid-19 pandemics remains constant [44]. Therefore, all the state variables are bounded during the Covid-19 pandemics.

Remark 3

In this research, it is assumed that vaccines generate complete immunity against coronavirus [45, 46]. Therefore, vaccinated cases move from the susceptible to the recovered compartment directly.

3 Estimation

In this section, the unknown parameters and noisy states of the system are estimated by the LSTM and EKF methods, respectively.

3.1 Parameters estimation by LSTM

In this section, the transmission rates of SIQHRE model, i.e., \(\Gamma =\{\alpha ,\gamma ,\lambda ,\theta ,\mu ,\phi \}\), are estimated by the LSTM method. In the LSTM method, the memorizing process is controlled by gates, and long-term dependency is prevented. As shown in Fig. 2, a typical LSTM unit has a cell state and three gates, including an input gate, a forget gate, and an output gate. These gates are some neural networks that can learn what information is relevant to keep or forget on the cell state during training. By means, LSTM can transfer associated information to the next unit of the long chain by these gates. Also, a gate is similar to a neural network layer, which contains different individual weights.

The cell state is the “memory” of the network and transfers relative information to the output of the unit by adding or removing information through forget and input gates as follows:

Fig. 2
figure 2

Schematic of the LSTM unit

The cell state is the “memory” of the network and transfers relative information to the output of the unit by adding or removing information through forget and input gates as follows:

$$\begin{aligned}&C_{t}=f_{t}* C_{t-1}+ i_{t}* \widetilde{C}_{t} \end{aligned}$$
(7)

where \(C_{t}, f_{t}, i_{t}\), and \(\widetilde{C}_{t}\) are current cell state value, current forget gate value, the output of sigmoid, and tangent hyperbolic functions in the input gate, respectively.

The input gate updates the cell state value. By means, previous hidden state and the current input pass into sigmoid and hyperbolic tangent functions as presented in Eqs. 8and 9.

$$\begin{aligned} i_{t}&={\mathrm {sigmoid}}(W_{i}.[h_{t-1},x_{t} ]+b_{i}) \end{aligned}$$
(8)
$$\begin{aligned} \widetilde{C}_{t}&=tanh (W_{c}.[h_{t-1},x_{t}]+b_{c}) \end{aligned}$$
(9)

where \(h_{t-1}\) and \(x_{t}\) present the previous hidden state and current input of the LSTM unit, respectively. The term \(W_{i}\) indicates the weighting matrix of the sigmoid operator between the input and output gates, and \(W_{C}\) is the weighting matrix of cell state. The terms \(b_{i}\), \(b_{c}\) are bias vectors related to the \(W_{i}\) , \(W_{C}\), respectively.

In the next step, the outputs are multiplied together to choose which information is more important to be kept. According to Eq. 10, the forget gate determines which data from the previous hidden state or the current input should be kept by a sigmoid function.

$$\begin{aligned}&f_{t}={\mathrm {sigmoid}}(W_{f}.[h_{t-1)},x_{t} ]+b_{f}) \end{aligned}$$
(10)

where \(W_{f}\) is the weighting matrix of the forget gate and \(b_{f}\) is the bias vectors related to \(W_{f}\). According to Eq. 11, in the output gate, the previous hidden state and the current input are passed into a sigmoid function and multiply with the \(tanh(C_{t})\) to decide which information should be transferred with the hidden state.

$$\begin{aligned}&h_{t}={\mathrm {sigmoid}}(W_{o}.[h_{t-1},x_{t} ]+b_{o})*tanh(C_{t}) \end{aligned}$$
(11)

where \(W_{o}\) is weighting matrix of output gate and \(b_{o}\) is the bias vector related to \(W_{o}\). The weighting matrix and bias vectors must be updated by gradient backpropagation to outperform. The parameter updated formulas are defined as follows:

$$\begin{aligned} \delta W&=\sum _{t=0}^{T}\delta {\mathrm {gate}}_{t}\bigotimes x_{t}+\sum _{t=0}^{T-1} \delta {\mathrm {gate}}_{t-1}\bigotimes h_{t} \end{aligned}$$
(12)
$$\begin{aligned} \delta b&=\sum _{t=0}^{T} \delta {\mathrm {gate}}_{t+1} \end{aligned}$$
(13)

where the variables \(gate_{t}\), W, and b are defined as follows:

$$\begin{aligned}&{\mathrm {gate}}_{t}=[C_{t} \quad i_{t} \quad f_{t} \quad O_{t} ]^{T} \end{aligned}$$
(14)
$$\begin{aligned}&W=[W_{c} \quad W_{i} \quad W_{f} \quad W_{o}]^{T} \end{aligned}$$
(15)
$$\begin{aligned}&b= [b_{c} \quad b_{i} \quad b_{f} \quad b_{o}]^{T} \end{aligned}$$
(16)

The estimated parameters vector is shown as, \({\widetilde{\Lambda }} =\{{\widetilde{\alpha }},{\widetilde{\gamma }},{\widetilde{\lambda }},{\widetilde{\theta }} ,{\widetilde{\mu }}\}\).

3.2 States estimation by EKF

In this section, unmeasurable states of the SIQHRE model are estimated by the EKF algorithm. There are two reasons why the EKF algorithm is necessary to be employed. First, accurate measurement of susceptible, quarantined, and healed groups is impossible. The second reason is that RT-PCR tests have the possibility of false-positive and false-negative results. Therefore, the measurement of infected cases is noisy and cannot be employed for real-time control purposes.

In order to employ the EKF algorithm, the mathematical description of the SIQHRE model is expressed as follows:

$$\begin{aligned}&x_{k+1}=f(x_{k},u_{k},w_{k })\quad {\mathrm {with}} \quad w(t)~N(0,Q) \end{aligned}$$
(17)
$$\begin{aligned}&y_{k}=h(x_{k},v_{k} ) \quad {\mathrm {with}} \quad v(t)~N(0,R) \end{aligned}$$
(18)

The terms \(w(t)\in \mathfrak {R}^{6\times 1}\) and \(v(t)\in \mathfrak {R}^{2\times 1}\) are process and measurement noises with the covariance matrixes Q and R, respectively. As a first step, it is necessary to compute the Jacobean matrices of the system as follows:

$$\begin{aligned}&A_{k}= \frac{\delta f({\hat{x}}_{k},u_{k},0)}{\delta x_{k} },H_{k}= \frac{ \delta h({\hat{x}}_{k+1},0)}{\delta x_{k}} \end{aligned}$$
(19)
Fig. 3
figure 3

EKF algorithm

According to Fig. 3, the EKF algorithm includes two phases, i.e., time update and measurement update. In the time-update phase, the predicted state estimate and covariance estimate are computed as follows:

$$\begin{aligned}&{\hat{x}}_(k+1)^{-}=f({\hat{x}}_{k},u_{k},0) \end{aligned}$$
(20)
$$\begin{aligned}&{\hat{P}}_{k+1}^{-}=A_{k} P_{k} A_{k}^{T}+Q_{k} \end{aligned}$$
(21)

In the measurement-update phase, the Kalman gain, updated-state estimate, and the updated covariance estimate are computed as follows:

$$\begin{aligned}&K= {\hat{P}}_{k+1}^{-} H_{k}^{T} (H_{k} {\hat{P}}_{k+1}^{-} H_{k}^{T}+R_{k} )^{-1} \end{aligned}$$
(22)
$$\begin{aligned}&{\hat{x}}_{k+1}={\hat{x}}_{k+1}^{-}+K(y_{k}-h({\hat{x}}_{k+1}^{-} )) \end{aligned}$$
(23)
$$\begin{aligned}&P_{k+1}=(I-KH_{k} ) {\hat{P}}_{k+1}^{-} \end{aligned}$$
(24)

It is important to note that the norm of EKF estimation error is bounded by a positive scalar, i.e., \(||{\hat{e}}_{\Omega }||\le \varepsilon _{\Omega }\), \(\Omega \in \{S,I,Q,R\}\), if disturbing noise and initial estimation are bounded. Therefore, the estimated variables are denoted as \({\hat{\Omega }}=\Omega _{n}+{\hat{e}}_{\Omega }\) [47], where \({\hat{\Omega }}\) and \(\Omega _{n}\) are estimated variables and nominal value of the estimated variables, respectively. In this research, the estimation error bounds (\(\varepsilon _{\Omega }\)) are obtained through trial and error.

Fig. 4
figure 4

Structure of closed-loop system

4 Control design

In this section, a sliding mode controller based on the LSTM and EKF estimators is developed to mitigate the spread of coronavirus in the Iranian and Russian societies. The structure of the closed-loop system is depicted in Fig. 4. In the next step, the dynamical equations of the SIQRHE system are rewritten based on LSTM and EKF estimators as follows:

$$\begin{aligned}&S(k+1)=S(k)- {\widetilde{\alpha }} (I_{n} (k)+{\hat{e}}_{I} (k))(S_{n} (k)+{\hat{e}}_{s} (k)) \end{aligned}$$
(25)
$$\begin{aligned}&I(k+1)= \ I(k)+ {\widetilde{\alpha }} (I_{n} (k)+{\hat{e}}_{I} (k))(S_{n} (k)+{\hat{e}}_{s} (k))\nonumber \\&\quad -({\widetilde{\gamma }}+{\widetilde{\lambda }})( I_{n}(k)\nonumber \\&\quad +{\hat{e}}_{I}(k))+u2 \end{aligned}$$
(26)
$$\begin{aligned}&Q(k+1)=Q(k)+{\widetilde{\lambda }}(I_{n} (k)+{\hat{e}}_{I} (k))-({\widetilde{\theta }}\nonumber \\&\quad +{\widetilde{\mu }})(Q_{n} (k)+{\hat{e}}_{Q} (k)) \end{aligned}$$
(27)
$$\begin{aligned}&H(k+1)= \ H(k)+{\widetilde{\gamma }} (I_{n} (k)+{\hat{e}}_{I} (k))\nonumber \\&\quad +{\widetilde{\theta }} (Q_{n} (k)+{\hat{e}}_{Q} (k))\nonumber \\&-({\widetilde{\phi }}+{\widetilde{\varepsilon }})H(k) \end{aligned}$$
(28)
$$\begin{aligned}&R(k+1)=R(k)+{\widetilde{\phi }} H(k)+{\widetilde{\mu }}(Q_{n} (k)+{\hat{e}}_{Q} (k))+u1 \end{aligned}$$
(29)
$$\begin{aligned}&E(k+1)=E(k)+{\widetilde{\varepsilon }} H(k)+u3 \end{aligned}$$
(30)

For the sake of simplicity, it is recommended to separate the nominal terms from uncertain terms. Therefore, the uncertain terms are written in the bracket as follows:

$$\begin{aligned}&S(k+1)= \ S(k)-{\widetilde{\alpha }} (I_{n} (k)+S_{n} (k))\nonumber \\&\quad -[{\widetilde{\alpha }}({\hat{e}}_{I} (k) {\hat{e}}_{s} (k)\nonumber \\&\quad +I_{n} (k){\hat{e}}_{s} (k)+S_{n} (k) {\hat{e}}_{I} (k))] \end{aligned}$$
(31)
$$\begin{aligned}&I(k+1)= \ I(k)+ {\widetilde{\alpha }} (I_{n} (k)+S_{n} (k))\nonumber \\&\quad +({\widetilde{\gamma }}+{\widetilde{\lambda }} ) I_{n} (k)\nonumber \\&\quad +[{\widetilde{\alpha }}({\hat{e}}_{I} (k){\hat{e}}_{s} (k)\nonumber \\&\quad \ +I_{n} (k){\hat{e}}_{s} (k) + S_{n} (k){\hat{e}}_{I} (k))\nonumber \\&\quad +({\widetilde{\gamma }}+{\widetilde{\lambda }} ){\hat{e}}_{I} (k)]+u2 \end{aligned}$$
(32)
$$\begin{aligned}&Q(k+1)=Q(k)+{\widetilde{\lambda }}(I_{n} (k))-({\widetilde{\theta }}\nonumber \\&\quad +{\widetilde{\mu }} )(Q_{n} (k))+[{\widetilde{\lambda }}{\hat{e}}_{I} (k)\nonumber \\&\quad -({\widetilde{\theta }}+{\widetilde{\mu }} ){\hat{e}}_{Q} (k)] \end{aligned}$$
(33)
$$\begin{aligned}&H(k+1)=H(k)+{\widetilde{\gamma }} (I_{n} (k))+{\widetilde{\theta }}(Q_{n} (k))\nonumber \\&\quad -({\widetilde{\phi }}+{\widetilde{\varepsilon }} )H(k)\nonumber \\&\quad +[{\widetilde{\gamma }}{\hat{e}}_{I} (k)+{\widetilde{\theta }}{\hat{e}}_{Q} (k)] \end{aligned}$$
(34)
$$\begin{aligned}&R(k+1)=R(k)+{\widetilde{\phi }} H(k)+{\widetilde{\mu }}(Q_{n} (k))\nonumber \\&\quad +[{\widetilde{\mu }}{\hat{e}}_{Q} (k)]+u1 \end{aligned}$$
(35)
$$\begin{aligned}&E(k+1)=E(k)+{\widetilde{\varepsilon }}H(k)+u3 \end{aligned}$$
(36)

4.1 Sliding mode control

The procedure of sliding mode controller design is initiated by defining the sliding manifolds as follows:

$$\begin{aligned}&S_{R} (k)= R(k)-R_{d} (k) \end{aligned}$$
(37)
$$\begin{aligned}&S_{I} (k)= I(k)-I_{d} (k) \end{aligned}$$
(38)
$$\begin{aligned}&S_{E} (k)= E(k)-E_{d} (k) \end{aligned}$$
(39)

where \(R_{d} (k),I_{d} (k)\) and \(E_{d} (k)\) are desired behavior of recovered, infected, and extinct compartments. The equivalent control laws are responsible for moving the state variables toward sliding manifolds, i.e., (\(S_{i} (k+1)-S_{i} (k)=0 for i\in \{R,I,E\}\)). To obtain the equivalent control laws, the forward difference of sliding manifolds is computed, and the nominal terms of Eqs. 3235 and 36 are substituted as follows:

$$\begin{aligned}&S_{R} (k+1) -S_{R} (k)\ = {\widetilde{\phi }} H(k)\nonumber \\&\quad + {\widetilde{\mu }} Q_{n} (k)+u_{R}-(R_{d} (k+1)-R_{d} (k))\nonumber \\&=0 \end{aligned}$$
(40)
$$\begin{aligned}&S_{I} (k+1)-S_{I} (k)= \ {\widetilde{\alpha }} (I_{n} (k) +S_{n} (k))\nonumber \\&\quad +({\widetilde{\gamma }}+{\widetilde{\lambda }} ) I_{n} (k)+u_{I}\nonumber \\&\quad -(R_{I} (k+1)-R_{I} (k))=0 \end{aligned}$$
(41)
$$\begin{aligned}&S_{E} (k+1)-S_{E} (k)={\widetilde{\varepsilon }} H(k)\nonumber \\&\quad +u_{E}-(R_{E} (k+1)-R_{E} (k))=0 \end{aligned}$$
(42)

After some mathematical calculations, the equivalent control laws are achieved as follows:

$$\begin{aligned}&u_{R}=-{\widetilde{\phi }} H(k)-{\widetilde{\mu }} Q_{n} (k)+(R_{d} (k+1)-R_{d} (k)) \end{aligned}$$
(43)
$$\begin{aligned}&u_{I}=-{\widetilde{\alpha }} (I_{n} (k)+S_{n} (k))-({\widetilde{\gamma }}\nonumber \\&\quad + {\widetilde{\lambda }} ) I_{n} (k)+(R_{I} (k+1)-R_{I} (k)) \end{aligned}$$
(44)
$$\begin{aligned}&u_{E}=-{\widetilde{\varepsilon }} H(k)+(R_{E} (k+1)-R_{E} (k)) \end{aligned}$$
(45)

In continuation, the state variables of the SIQHRE model move toward origin along the sliding manifolds by applying switching control laws as follows:

$$\begin{aligned}&u1=u_{R}- K_{R}(k) sgn ( S_{R} (k)) \end{aligned}$$
(46)
$$\begin{aligned}&u2=u_{I}- K_{I}(k) sgn ( S_{I} (k)) \end{aligned}$$
(47)
$$\begin{aligned}&u3=u_{E}- K_{E}(k) sgn ( S_{E} (k)) \end{aligned}$$
(48)

where \(K_{R}(k),K_{I}(k)\), and \(K_{E}(k)\) are positive time-dependent switching gains. The stability of the closed-loop system is guaranteed through the proper selection of switching gains. For this purpose, a rigorous stability analysis is proposed in the next section.

4.2 Stability analysis

The positive Lyapunov candidate function is considered as follows:

$$\begin{aligned}&V_{R} (k)=\frac{1}{2} (S_{R}^{2} (k)+S_{I}^{2} (k)+S_{E}^{2} (k)) \end{aligned}$$
(49)

Taking the forward difference of both sides of the above equation and using reachability condition yield:

$$\begin{aligned} \begin{aligned}&V_{R} (k+1)-V_{R} (k)={} \ S_{R} (k)(S_{R} (k+1)\\&\quad -S_{R} (k))+S_{I} (k)(S_{I} (k+1)\\&\quad \ -S_{I} (k) +S_{E} (k)(S_{E} (k+1)-S_{E} (k))\le -\zeta _{R} \\&\quad \ |S_{R} (k)| -\zeta _{I} |S_{R} (k)| -\zeta _{E} |S_{E} (k)| \end{aligned} \end{aligned}$$
(50)

where \(\zeta _{i}\ge 0,i\in {R,I,E}\) are design parameters. Substitution of Eqs. 32, 35 and 36 into the above equation yields:

$$\begin{aligned} \begin{aligned}&V_{R} (k+1)-V_{R} (k)={} \ S_{R} (k)({\widetilde{\phi }}H(k)\\&\quad +{\widetilde{\mu }}Q_{n} (k)+{\widetilde{\mu }} {\hat{e}}_{Q} (k) \\&\quad \ +u1-(R_{d} (k+1)-R_{d} (k)))+S_{I} (k)({\widetilde{\alpha }} (I_{n} (k)\\&\quad \ +S_{n} (k))+({\widetilde{\gamma }}\\&\quad +{\widetilde{\lambda }} ) I_{n} (k)+u_{I}-(R_{I} (k+1)\\&\quad \ -R_{I} (k)))+S_{E} (k)({\widetilde{\varepsilon }}H(k)+u_{E}-(R_{E} (k+1)\\&\quad \ -R_{E} (k)))\le -\zeta _{R} |S_{R} (k)|-\zeta _{I} |S_{R} (k)| \\&\quad \ -\zeta _{E} |S_{E} (k)| \end{aligned} \end{aligned}$$
(51)

The substitution of Eqs. 43-45 into Eq. 48 yields:

$$\begin{aligned} \begin{aligned}&V_{R} (k+1)-V_{R} (k)={} \ S_{R} (k)({\widetilde{\mu }} \\&\quad {\hat{e}}_{Q} (k) -K_{R}(k) sgn ( S_{R} (k)))+S_{E} (k)\\&\quad \ ({\widetilde{\alpha }}({\hat{e}}_{I} (k){\hat{e}}_{s} (k)\\&\quad +I_{n} (k) {\hat{e}}_{s} (k) + S_{n} (k) {\hat{e}}_{I} (k))\\&\quad \ +({\widetilde{\gamma }}+{\widetilde{\lambda }} ) \\&\quad {\hat{e}}_{I} (k)- K_{I}(k) sgn ( S_{I} (k))+S_{I} (k)\\&\quad \ (- K_{E}(k) sgn ( S_{E} (k))\le -\zeta _{R}|S_{R} (k)|\\&\quad \ -\zeta _{I} |S_{R} (k)| -\zeta _{E} |S_{E} (k)| \end{aligned} \end{aligned}$$
(52)

After some mathematical simplification, the above equation can be written as follows:

$$\begin{aligned} \begin{aligned}&{} \ V_{R} (k+1)-V_{R} (k)\le {\widetilde{\mu }}\varepsilon _{Q} |S_{R} (k)| \\&\quad -K_{R}(k) |S_{R} (k)|+|S_{E} (k)| (-{\widetilde{\alpha }}\varepsilon _{I} \varepsilon _{S}\\&\ (+({\widetilde{\gamma }}+{\widetilde{\lambda }}) \varepsilon _{I}\\&\quad +{\widetilde{\alpha }}\varepsilon _{S} |I_{n} (k)|\\&\quad +\varepsilon _{I} |S_{n} (k)|))- K_{I}(k) |S_{I} (k)||\\&\ - K_{E}(k) |S_E (k) \le -\zeta _{R} |S_{R} (k)|\\&\quad -\zeta _{I} |S_{R} (k)| -\zeta _{E} |S_{E} (k)| \end{aligned} \end{aligned}$$
(53)

According to the above equations, the state variables of the SIQHRE model converge exponentially to the desired values, and also, Lyapunov function is negative definite if switching gains are selected as follows:

$$\begin{aligned}&K_{R}(k)\ge \zeta _R+{\widetilde{\mu }}\varepsilon _{Q} \end{aligned}$$
(54)
$$\begin{aligned}&K_{I}(k)\ge \zeta _{I}+-{\widetilde{\alpha }}\varepsilon _{I} \varepsilon _{S}\nonumber \\&\quad +({\widetilde{\gamma }}+{\widetilde{\lambda }} ) \varepsilon _{I}\nonumber \\&\quad +{\widetilde{\alpha }}(\varepsilon _{S} |I_{n} (k)|+\varepsilon _{I} |S_{n} (k)|) \end{aligned}$$
(55)
$$\begin{aligned}&K_{E}(k)\ge \zeta _{E} \end{aligned}$$
(56)

It should be mentioned that proper selection of time-varying switching gains reduces the effect of system uncertainties on the obtained results in the next section.

Table 1 List of parameters

5 Simulation results

In this section, the performances of the proposed algorithms are analyzed. The parameters of the SIQHRE model are listed in Table 1. In this study, the official confirmed data provided by the Iranian and Russian ministries of health from 1 August 2020 to 20 March 2021 are employed [48,49,50]. In the proposed LSTM algorithm, data from 1 August 2020 to 9 January 2021 are utilized for training and validation purposes. The data from 10 January 2021 to 20 March 2021 are used for test purposes.

Fig. 5
figure 5

Estimation of recovery rate by LSTM method in Iran

Fig. 6
figure 6

Estimation of recovery rate by LSTM method in Russia

Fig. 7
figure 7

Estimation of mortality rate by the LSTM method in Iran

Fig. 8
figure 8

Estimation of mortality rate by the LSTM method in Russia

In Figs. 5, 6, the recovery rate is estimated by the LSTM method. Clearly, the results of the proposed LSTM method approximately fit with the actual data. It should be mentioned that during this period, the recovery rate in Iran experiences an incremental slope. In this estimation, the root-mean-square error (RMSE) and Pearson correlation are 0.1479 and 0.9904, respectively.

The actual and estimated mortality rates are compared in Figs. 7, 8. According to these figures, the mortality rate experiences a decreasing slope in Iran. During this period, the mortality rate shows an increasing slope in Russia. The proposed LSTM algorithm successfully tracks the actual values. The RMSE and Pearson correlation are 0.1394 and 0.9964, respectively.

Fig. 9
figure 9

Estimation of infection rate by LSTM method in Iran

Fig. 10
figure 10

Estimation of infection rate by LSTM method in Russia

In Figs. 9, 10, the actual and estimated data for the rate of infected cases are compared. Despite sharp variation in the actual data of Iran, the LSTM algorithm successfully estimates the actual data. In this period, the infection rate experiences a decreasing slope in Russia. The RMSE and Pearson correlation are 0.05445 and 0.9922, respectively.

Fig. 11
figure 11

Estimation of hospitalization rate by LSTM method in Iran

Fig. 12
figure 12

Estimation of hospitalization rate by LSTM method in Russia

In Figs. 11, 12, the estimated rate of hospitalization is compared with the actual data. According to these figures, in spite of sharp variations in the actual data of Iranian and Russian countries, the LSTM algorithm is able to estimate the actual data with high precision. The RMSE and Pearson correlation indexes are 0.2753 and 0.9490, respectively.

In this section, the performance of the proposed EKF algorithm is investigated. As mentioned previously, the measurable state variables are hospitalized cases (H), who received medical service in the governor hospitals, and confirmed extinct cases (E). The rest of the state variables are contaminated with white Gaussian noises. In this study, the covariance of measurement and process noises and also the initial covariance estimation are defined as \(R=0.01 I_{2}, Q= I_{6},\) and \(P_{0}^{+}=10I_{6}\), respectively. In Figs. 13, 14, the tracking errors of suspected, infected, quarantined, and recovered cases are depicted. Clearly, the immeasurable state variables converge successfully to the actual data under the proposed EKF algorithm

Fig. 13
figure 13

Estimation errors of the EKF algorithm in Iran

Fig. 14
figure 14

Estimation errors of the EKF algorithm in Russia

Finally, the performance of the proposed control algorithm is studied. At first, the desired behavior of the system’s outputs, i.e., infected, extinct, and recovered, is defined as follows:

$$\begin{aligned}&I_{d} (k)= (I_{0}-I_{f} )\frac{cos(w_{I} k)}{k+1}+I_{f} \end{aligned}$$
(57)
$$\begin{aligned}&E_{d}(k)= (E_{0} - E_{f}) {\mathrm {exp}}(-w_{E} k)+E_{f} \end{aligned}$$
(58)
$$\begin{aligned}&R_{d}(k)= (R_{0} - R_{f}) {\mathrm {exp}}(w_{R} k)+R_{f} \end{aligned}$$
(59)

where \(I_{f}, E_{f}\), and \(R_{f}\) are desired final number of infected, extinct, and recovered compartments, respectively. As presented in Table 1, the parameters \(I_{0}, E_{0}\), and \(R_{0}\) are the initial number of these compartments and the parameters \(w_{I}, w_{E}\), and \(w_{R}\) are oscillatory rates and take positive values.

Fig. 15
figure 15

The effect of medical treatment on the daily deaths in Iran

Fig. 16
figure 16

The effect of medical treatment on the daily deaths in Russia

According to Figs. 15, 16, Covid-19 daily deaths in Iran and Russia experience a decreasing slope under improvement of medical treatment. In other words, positive actions in the medical service system such as increasing hospital capacity, improving hospitals equipment, and employing fresh nurses and doctors lead to decreasing mortality in a short time. The proposed controller can reduce the Covid-19 daily death to less than ten people in less than one month in Iran and Russia.

Fig. 17
figure 17

The effects of social distance and facial mask on daily infected cases

Fig. 18
figure 18

The effects of social distance and facial mask on daily infected cases

In Figs. 17, 18, the daily confirmed infected cases under social distance and facial masks policies are presented. According to these figures, daily Covid-19-infected cases become less than 500 people in less than one month under the proposed controller. However, after the peak of infected cases, a steep decrease is observed in daily confirmed infected cases. It is concluded that social distance and facial masks are effective policies to reduce the daily infection rate. Since daily infection rate has a direct relationship with the mortality rate of Covid-19, social distance and facial mask result in the reduction of Covid-19 deaths in the long term.

Fig. 19
figure 19

The effect of vaccination on the immunity of Iranian society against Covid-19 pandemic

Fig. 20
figure 20

The effect of vaccination on the immunity of Iranian society against Covid-19 pandemic

In the last simulation, the effect of vaccination on the immunity of people against novel Covid-19 is investigated. During this simulation, it is assumed that vaccination was initiated on 10 January 2021, and also, vaccinated people have complete immunity against Covid-19. According to Figs. 19, 20, from 10 January 2021 to 25 February 2021, the vaccination has no tangible effect on the Iranian and Russian societies. Between 25 February 2021 and 20 March 2021, vaccination produces herd immunity for a number of 80 million people in Iran. Due to the higher population of the Russian community, the herd immunity occurred later than Iranian community. It is concluded that vaccination as an effective action against the Covid-19 pandemic is able to produce herd immunity in the long term. However, it should be mentioned that social distance and facial masks are obligatory during global vaccination.

6 Conclusion

In this paper, a comprehensive solution is proposed to mitigate the spread of Covid-19 in the Iranian and Russian societies. At first, the transmission flow of Covid-19 is expressed by the SIQHRE model. The transmission rates of the system are estimated by the LSTM algorithm. Then, the EKF algorithm is employed to estimated state variables of the system in the presence of measurement and process noises. In continuation, a robust sliding mode controller is developed to control the spread of Covid-19 under vaccination, social distance and facial masks, and medical treatment. The stability of the closed-loop system is analyzed through Lyapunov theorems. This research has the potential to be extended in different ways. For example, it is recommended to consider the possibility of infection after vaccination. The results of this research help policy-makers to adopt the best decisions in a global fight against Covid-19. However, it should be mentioned that there are several obstacles to eliminating the Covid-19 pandemic, including the emergence of novel variants of Covid-19, such as SARS-CoV-2, and also different levels of vaccination protection against new variants of Covid-19.