Model Predictive Control of COVID-19 Pandemic with Social Isolation and Vaccination Policies in Thailand

Abstract

This study concerns the COVID-19 pandemic in Thailand related to social isolation and vaccination policies. The behavior of disease spread is described by an epidemic model via a system of ordinary differential equations. The invariant region and equilibrium point of the model, as well as the basic reproduction number, are also examined. Moreover, the model is fitted to real data for the second wave and the third wave of the pandemic in Thailand by a sum square error method in order to forecast the future spread of infectious diseases at each time. Furthermore, the model predictive control technique with quadratic programming is used to investigate the schedule of preventive measures over a time horizon. As a result, firstly, the plan results are proposed to solve the limitation of ICU capacity and increase the survival rate of patients. Secondly, the plan to control the outbreak without vaccination shows a strict policy that is difficult to do practically. Finally, the vaccination plan significantly prevents disease transmission, since the populations who get the vaccination have immunity against the virus. Moreover, the outbreak is controlled in 28 weeks. The results of a measurement strategy for preventing the disease are examined and compared with a control and without a control. Thus, the schedule over a time horizon can be suitably used for controlling.

1. Introduction

COVID-19 is a new disease that appeared at the end of 2019 and is caused by the virus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The virus spreads to others if they come into contact with respiratory droplets or aerosols of infected people through coughing, sneezing, or talking, as well as through close contact or sharing items with an infected person [1,2]. Infected individuals can present symptoms or be asymptomatic [3]. Some studies have suggested the proportion of infected people who show symptoms is eighty percent, with asymptomatic individuals accounting for twenty percent [4]. Fever, dry cough, and fatigue are among the most common symptoms of people who are infected and show symptoms [5]. After people get the virus, they can detect the infection and show symptoms after 2–14 days, which is the incubation period of COVID-19 [6].

The COVID-19 pandemic began in December 2019 in Wuhan, China. It continued to spread around the world and has killed millions of people, as shown in the report of the World Health Organization (WHO) [7]. The total number of infected people and deaths worldwide was 210,805,121 and 4,416,421, respectively. The data were last updated on 19 August 2021. Aside from the epidemic having a negative impact on the health of people, it has affected education, the economy, human life, and others. For example, laying off employees because of lockdown measures that forced businesses to close, testing the effectiveness of online education during the pandemic, and food security. Therefore, locking down, social distancing, wearing a face mask, and working from home were preventive measures to be used to prevent the transmission of the COVID-19 pandemic in many countries. For example, a survey of more than 2000 adults in the United States reported that 72 percent intend to continue wearing face masks in public, 88 percent intend to avoid crowds, and 50 percent intend to frequently wash their hands and use hand sanitizer [8]. A vaccine to prevent COVID-19 infection was developed, and the mRNA vaccine, the Pfizer/BioNTech vaccine, was the first to receive emergency validation from the World Health Organization (WHO). Now, the WHO has evaluated the following vaccines: AstraZeneca/Oxford, Johnson and Johnson, Moderna, Pfizer/BioNTech, Sinopharm, and Sinovac, against COVID-19 infection. These vaccines were used in many countries to reduce the number of infected people [9,10].

The COVID-19 epidemic has already spread to Thailand in three waves. The first wave of the outbreak began in January 2020, with a peak of 188 infected patients per day, and was finally brought under control in March, while the second wave started in November 2020, with a rise in infected populations that reached 800 cases a day, before decreasing in March 2021. Furthermore, the third wave started in April 2021. On 19 August 2021, the number of newly infected cases and death cases each day was higher than in the preceding two waves, with 20,902 cases and 235 cases, respectively. The data were reported by Thailand’s Ministry of Public Health’s Department of Disease Control and the World Health Organization (WHO) [11,12]. The total number of infected cases in Thailand was 989,859, with 8586 deaths. Social isolation policies, such as social distancing and the use of face masks, were used to control the epidemic in the three waves. Furthermore, a vaccine against COVID-19 infection became critical for the control of disease spread in the third wave, which was mostly limited to healthcare workers and relied primarily on the Sinovac Biotech vaccine. The vaccine was imported from China. The advantage of this vaccine is the reduction in the risk of death from infection. In contrast, the efficiency of preventing the virus is lowest compared with other kinds of vaccine, such as the Oxford–AstraZeneca vaccine. The Oxford–AstraZeneca vaccine was planned to be administered to the majority of the population to prevent COVID-19 infections [13]. Thus, this vaccine is investigated in this work.

In order for people to live a normal life, the outbreak should be resolved. Furthermore, efficient management can mitigate damage during a pandemic. Therefore, the application of mathematical fields is used to study the pandemic and compute a suitable plan for preventive measures.

Mathematical modeling is one of the applications that can be used for solving real-world problems. It can be applied to study and predict the behavior of outbreaks via a system of nonlinear differential equations, such as the models of influenza, Ebola, and coronavirus diseases by [4,14,15]. Now, epidemic models have been developed. Chancharoenthana Wiwat et al. used the $SIR$ model (Susceptible–Infected–Recovered) to study the COVID-19 pandemic and examined the contacts and the recovery of exposed and infected individuals [16]. After that, the $SIR$ model was extended and then became the $SEIR$ model (Susceptible–Exposed–Infected–Recovered). Natcha C Jitsuk et al. examined the impact of Thailand’s Songkran festival on COVID-19 transmission by using this model [17]. Thus, the incubation period of the disease was considered in the model. Moreover, the kinds of populations were determined and then the model was improved. The contact between and recovery of exposed and infected populations, symptomatic and asymptomatic infected populations, were focused on by T D Frank et al. and Robert C. Reiner Jr et al. [18,19]. In addition, Saulo B. Bastos et al. studied the number of death cases and predicted the early evolution of the COVID-19 pandemic in Brazil [20]. Furthermore, a quarantine was applied in an epidemic model to establish its association between infection and the control of disease spread. Enahoro Iboi et al. organized the classical model to assess the effectiveness of public health education campaigns during the COVID-19 pandemic in the United States. This model concerns the isolation and non-isolation of people, as well as the number of hospitalized and non-hospitalized individuals [21], and Pakwan Riyapan et al. concentrated on the isolation of infected people. Calistus N. Ngonghala et al. focused on the isolation of non-infected and infected people, and looked into using the quarantine model to resolve the COVID-19 pandemic. They used sensitivity analysis to investigate the interaction of the outbreak and the parameters of transmission, quarantine, and contact tracing [4,22].

The $SEIR$ model was depicted by fractional models in order to investigate disease transmission. For example, Saheed O. Akindeinde et al. used the fractional model to study the interaction of humans and pathogens, Olaniyi Iyiola et al. used the model to focus on the quarante of exposed and infected populations [23,24], and Isaac Owusu-Mensah et al. organized the model to investigate the rates of transmission, testing, and the transition of asymptomatic to symptomatic infected people and then continued to examine the impact of the rates on the outbreak [25]. Moreover, social isolation measures were studied in the model of the COVID-19 pandemic in Wuhan, Italy, and the United States, as well as the first and second waves of the outbreak in Thailand [4,26,27,28,29].

The optimum control theory is a method used to solve dynamic problems. It is a method used for the minimization of the cost function while computing control variables, for instance, applying optimal control of continuous-time and discrete-time in epidemic models. Pontryagin’s maximum was used to study the swine influenza disease, and model predictive control (MPC) was used to solve the Atopic Dermatitis (AD) model [14,30]. Moreover, quadratic programming is also being studied for use in the solving of constraint problems [31,32].

In this work, the mathematical model was developed from the model by Calistus N.Ngonghala et al. in order to study the second and third waves of the pandemic in Thailand [4]. The model focused on social distancing and vaccination policies, as well as the parameters of transmission, quarantine, and contact tracing. Model predictive control was used to solve discrete-time models. The technique provided the results of controlled variables over a time horizon. The results are suitable for a large group of people to use in reality when compared with the results of the continuous-time model. The idea of the technique was to compute control variables that can control the output of the system to reach a set point.

The research constructed a mathematical model to examine the COVID-19 pandemic concerning prevention strategies. The model, as well as the model’s analysis, invariant regions, equilibrium point, and the basic reproduction number, are all presented in Section 2. Section 3 investigates model predictive control, including linearization and discretization algorithms. In addition, Section 4 exhibits the curve fitting of the COVID-19 pandemic in Thailand and the pandemic’s prevention strategies. In the discussion and conclusion sections, the work is summarized.

2. Mathematical Model

The mathematical model of the COVID-19 pandemic is described in this section, followed by an analysis of the model that includes the invariant region, equilibrium points, and the basic reproduction number.

The disease transmission of the COVID-19 pandemic is illustrated by the flow diagram in Figure 1. The nine components in the diagram represent the kinds of populations. Non-quarantined susceptible ($S_u$), quarantined susceptible ($S_{sq}$), who are quarantined in state quarantine, non-quarantined exposed ($E_u$), quarantined exposed ($E_q$), exposed individuals who are infected but do not show symptoms and can not transmit infection, symptomatically-infectious ($I_s$), who are infected and show the symptoms of COVID-19, asymptomatically-infectious ($I_a$), who infected but exhibiting mild or no clinical symptoms of COVID-19, isolated infectious ($I_i$), who isolated in hospitals and state quarantine and isolated infections require intensive care unit ($I_{icu}$), as well as recovered individuals (R) at any time t. Moreover, the total population (N) is defined by

$$N\left(t\right)=S_u\left(t\right)+S_{sq}\left(t\right)+E_u\left(t\right)+E_q\left(t\right)+I_s\left(t\right)+I_i\left(t\right)+I_a\left(t\right)+I_{icu}\left(t\right)+R\left(t\right).$$

Moreover, the model is given as the system of a nonlinear differential equation in the following:

$$\begin{array}{cc}\hfill {\dot{S}}_u& =\mathsf{\Lambda }-((1-p)+(1-q)p+qp)\lambda S_u\left(t\right)+{\psi }_q{S}_{sq}\left(t\right)-\mu S_u\left(t\right)-{\epsilon }_v{u}_v{S}_u\left(t\right),\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill {\dot{S}}_{sq}& =(1-p)\lambda S_u\left(t\right)-({\psi }_q+\mu )S_{sq}\left(t\right)-{\epsilon }_v{u}_v{S}_{sq}\left(t\right),\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill {\dot{E}}_u& =(1-q)p\lambda S_u\left(t\right)-({\sigma }_u+{\alpha }_u+\mu )E_u\left(t\right),\hfill \end{array}$$

(1c)

$$\begin{array}{cc}\hfill {\dot{E}}_q& =qp\lambda S_u\left(t\right)+{\alpha }_u{E}_u\left(t\right)-({\sigma }_q+\mu )E_q\left(t\right),\hfill \end{array}$$

(1d)

$$\begin{array}{cc}\hfill {\dot{I}}_s& ={\delta }_u{\sigma }_u{E}_u\left(t\right)-(r_s+{\varphi }_s+m_s+\mu )I_s\left(t\right),\hfill \end{array}$$

(1e)

$$\begin{array}{cc}\hfill {\dot{I}}_i& =\gamma {\sigma }_q{E}_q\left(t\right)+{\varphi }_s{I}_s\left(t\right)+{\sigma }_a{I}_a\left(t\right)-(r_i+v_i+m_i+\mu )I_i\left(t\right),\hfill \end{array}$$

(1f)

$$\begin{array}{cc}\hfill {\dot{I}}_a& =(1-{\delta }_u){\sigma }_u{E}_u\left(t\right)+(1-\gamma ){\sigma }_q{E}_q\left(t\right)-(r_a+{\sigma }_a+m_a+\mu )I_a\left(t\right),\hfill \end{array}$$

(1g)

$$\begin{array}{cc}\hfill {\dot{I}}_{icu}& =v_i{I}_i\left(t\right)-(r_{icu}+m_{icu}+\mu )I_{icu}\left(t\right),\hfill \end{array}$$

(1h)

$$\begin{array}{cc}\hfill \dot{R}& =r_s{I}_s\left(t\right)+r_i{I}_i\left(t\right)+r_a{I}_a\left(t\right)+r_{icu}{I}_{icu}\left(t\right)-\mu R\left(t\right)+{\epsilon }_v{u}_v{S}_u\left(t\right)+{\epsilon }_v{u}_v{S}_{sq}\left(t\right),\hfill \end{array}$$

(1i)

where force of infection, $\lambda$, and initial condition, $IC$, are defined by

$$\lambda =\frac{{\beta }_{so}(1-u_{so})(1-u_M)(I_s\left(t\right)+{\eta }_a{I}_a\left(t\right))}{S_u\left(t\right)+E_u\left(t\right)+I_s\left(t\right)+I_a\left(t\right)+R\left(t\right)},$$

(2)

and

$$IC={\left[S_u\left(0\right)S_{sq}\left(0\right)E_u\left(0\right)E_q\left(0\right)I_s\left(0\right)I_i\left(0\right)I_a\left(0\right)I_{icu}\left(0\right)R\left(0\right)\right]}^T.$$

(3)

In addition, the number of deaths in the population (D) is

$$\dot{D}=m_s{I}_s\left(t\right)+m_i{I}_i\left(t\right)+m_a{I}_a\left(t\right)+m_{icu}{I}_{icu}\left(t\right).$$

(4)

The parameters $u_{so}$, $u_M$ and $u_v$ represent the effect of social distancing, the use of face masks and vaccination, $(0\le u_{so},u_M,u_v\le 1)$, which are set as constant, while $u_v$ is set as zero. The parameter $\mathsf{\Lambda }$ is the recruitment rate, $\mu$ is the natural mortality rate and $p\lambda$ is the rate of getting the virus with the probability of infection, p. In the susceptible classes, the non-quarantined populations are quarantined at a rate $(1-p)\lambda$, and they are released from quarantine and reversed to the non-quarantine class at a rate ${\psi }_q$. Moreover, they move to the exposed classes with and without quarantine at the rates $qp\lambda$ and $(1-q)p\lambda$, respectively, with the rate of quarantined, q. In the exposed class, there is no disease transmission because it is during the incubation stage of the virus. The incubation periods of non-quarantine and quarantine are defined by ${\sigma }_u$ and ${\sigma }_q$, respectively. The non-quarantined populations are quarantined at the rate of ${\alpha }_u$, contact tracking rate, and they move to the symptomatic and asymptomatic infectious classes at the rates ${\delta }_u{\sigma }_u$ and $(1-{\delta }_u){\sigma }_u$, respectively, at the end of the incubation period. In addition, the quarantined populations also move to the symptomatic and asymptomatic infectious classes at the rates $\gamma {\sigma }_q$ and $(1-\gamma ){\sigma }_q$, respectively, at the end of the incubation period. In infected classes, the symptomatic and asymptomatic infectious populations are isolated by contact tracking at the rate ${\varphi }_s$ and ${\sigma }_a$, respectively. Furthermore, the isolated people move to the intensive care unit classes at the rate $v_i$. Next, the recovery rates and mortality rates due to the COVID-19 pandemic of the infected classes, $I_s,I_i,I_a$, and $I_{icu}$, are defined by $r_s,r_i,r_a,r_{icu}$ and $m_s,m_i,m_a,m_{icu}$. Moreover, the susceptible population will have immunity against the disease at the rate ${\epsilon }_v{u}_v$ with vaccine efficacy, ${\epsilon }_v$.

2.1. Invariant Region

Invariant region is the method used for the study of the non-negative of all components of the model (1); it is shown in the following Theorem.

Theorem 1.

Suppose that the initial value $IC$ in Equation (3) is all non-negative. Then, the solution of model (1) that start with these initial values will remain non-negative for all time $t>0$. Furthermore, all solutions of the model (1) are bounded.

Proof of Theorem 1.

Suppose $\delta \in min\{m_s,m_i,m_a,m_{icu}\}$ and $d\in min\{\mu ,\delta \}$.

$$\begin{array}{cc}\hfill \dot{N}& =\mathsf{\Lambda }-\mu N-m_s{I}_s-m_i{I}_i-m_a{I}_a-m_{icu}{I}_{icu}\hfill \\ \hfill & \le \mathsf{\Lambda }-\mu (S_u+S_{sq}+E_u+E_q+R)-(\mu +\delta )(I_s+I_i+I_a+I_{icu})\hfill \\ \hfill & \le \mathsf{\Lambda }-d(S_u+S_{sq}+E_u+E_q+R)-d(I_s+I_i+I_a+I_{icu})\hfill \\ \hfill & \le \mathsf{\Lambda }-dN,\hfill \end{array}$$

then $lim{sup}_{t\to \infty }N\left(t\right)\le \frac{\mathsf{\Lambda }}{d}$, namely all solutions are non-negative. □

Theorem 2.

The region $\mathsf{\Gamma }=\{(S_u\left(t\right),S_{sq}\left(t\right),E_u\left(t\right),E_q\left(t\right),I_s\left(t\right),I_i\left(t\right),I_a\left(t\right),I_{icu}\left(t\right),R\left(t\right))\in {\mathbb{R}}_{+}^9:N\left(t\right)\le \frac{\mathsf{\Lambda }}{d}\}$ is positively invariant for the model (1) for every non-negative initial condition in ${\mathbb{R}}_{+}^9$.

From the above, the components of the model in Equation (1) satisfy Theorems 1 and 2; then, they are non-negative.

2.2. The Equilibrium Points

The equilibrium points of the model (1) yield two points: a disease-free equilibrium point ($E_{dfe}$) and an endemic equilibrium point ($E_{ee}$), which are given by the following equations:

• Disease-free equilibrium point $E_{dfe}$:

  $$\begin{array}{c}\hfill E_{dfe}=\left(S_u,S_{sq},E_u,E_q,I_s,I_i,I_a,I_{icu},R\right)=\left(\frac{\mathsf{\Lambda }}{\mu },0,0,0,0,0,0,0,0\right).\end{array}$$

  (5)
• Endemic equilibrium point $E_{ee}$:

  $$\begin{array}{c}\hfill E_{ee}=\left(S_u^{*},S_{sq}^{*},E_u^{*},E_q^{*},I_s^{*},I_i^{*},I_a^{*},I_{icu}^{*},R^{*}\right),\end{array}$$

  (6)

  where

  $$\begin{array}{cc}\hfill S_u^{*}& =\frac{\mathsf{\Lambda }+{\psi }_q{S}_{sq}^{*}}{A_1{\lambda }^{*}+\mu },S_{sq}^{*}=\frac{(1-p){\lambda }^{*}{S}_u^{*}}{A_2},E_u^{*}=\frac{(1-q)p{\lambda }^{*}{S}_u^{*}}{A_3},E_q^{*}=\frac{qp{\lambda }^{*}{S}_u^{*}+{\alpha }_u{E}_u^{*}}{A_4},\hfill \\ \hfill I_s^{*}& =\frac{{\delta }_u{\sigma }_u{E}_u^{*}}{A_5},I_i^{*}=\frac{\left(\gamma {\sigma }_q{E}_q^{*}+{\varphi }_s{I}_s^{*}+{\sigma }_a{I}_a^{*}\right)}{A_6},I_a^{*}=\frac{A_7{E}_u^{*}+E_q^{*}{\sigma }_q(1-\gamma )}{A_8},I_{icu}^{*}=\frac{v_i{I}_i^{*}}{A_9},\hfill \\ \hfill R^{*}& =\frac{\left(r_s{I}_s^{*}+r_i{I}_i^{*}+r_a{I}_a^{*}+r_{icu}{I}_{icu}^{*}\right)}{\mu },{\lambda }^{*}=\frac{A_0(I_s^{*}+{\eta }_a{I}_a^{*})}{S_u^{*}+E_u^{*}+I_s^{*}+I_a^{*}+R^{*}},\hfill \end{array}$$

  and $$\begin{gathered} A_0={\beta }_{so}(1-u_{so})(1-u_M),A_1=(1-p)+(1-q)p+qp,A_2={\psi }_q+\mu , \\ {A}_3={\sigma }_u+{\alpha }_u+\mu ,A_4={\sigma }_q+\mu ,A_5=r_s+{\varphi }_s+m_s+\mu ,A_6=r_i+v_i+m_i+\mu , \\ {A}_7=(1-{\delta }_u){\sigma }_u,A_8=r_a+{\sigma }_a+m_a+\mu ,A_9=r_{icu}+m_{icu}+\mu . \end{gathered}$$

2.3. The Basic Reproduction Number ${\mathcal{R}}_0$

The basic reproduction number is used to explain disease spread, the effectiveness of contact rate, the risk of infection, and the length of infectiousness by comparing the average number of newly infected cases with the number of infected cases in a community. The basic reproduction number of a model can be investigated by using the next-generation matrix method.

The method focuses on using the disease-free equilibrium point and the Jacobian matrix to apply a spectral radius theory to the infected components of a model. So, a subsystem of the model (1) concerning the states-at-infection form is the models (1c)–(1h). It can be rewritten as

$$\begin{array}{c}\hfill \dot{\mathbf{x}}=F\left(\mathbf{x}\right)-V\left(\mathbf{x}\right),\end{array}$$

(7)

where

$$\begin{array}{c}\hfill F\left(\mathbf{x}\right)=\left[\begin{array}{c}-\frac{A_0{S}_{u}p\left(I_a{\eta }_a+I_s\right)\left(q-1\right)}{E_u+I_a+I_s+R+S_u}\\ \frac{A_0{S}_{u}pq\left(I_a{\eta }_a+I_s\right)}{E_u+I_a+I_s+R+S_u}\\ 0\\ 0\\ 0\\ 0\end{array}\right],V\left(\mathbf{x}\right)=\left[\begin{array}{c}{A}_3{E}_u\\ A_4{E}_q-E_u{\alpha }_u\\ -E_u{\delta }_u{\sigma }_u+A_5{I}_s\\ -\gamma E_q{\sigma }_q+A_6{I}_i-I_a{\sigma }_a-I_s{\varphi }_s\\ -A_7{E}_u+A_8{I}_a+E_q{\sigma }_q\left(\gamma -1\right)\\ A_9{I}_{icu}-I_i{v}_i\end{array}\right],\end{array}$$

and $\mathbf{x}={\left[E_u{E}_q{I}_s{I}_i{I}_a{I}_{icu}\right]}^T$. Next, the spectral radius of the matrix, $J_F\left(E_{dfe}\right)J_V^{-1}\left(E_{dfe}\right)$, is defined by

$$\begin{array}{c}\hfill \rho \left(J_F\left(E_{dfe}\right)J_V^{-1}\left(E_{dfe}\right)\right)=\underset{i=1,2,\cdots ,n_e}{max}\left|e_i\right|,\end{array}$$

(8)

where $e_j$ are eigenvalues of the matrix and $n_e$ is number of the eigenvalues [33,34]. The matrix $J_F\left(E_{dfe}\right)$ and $J_V\left(E_{dfe}\right)$ are Jacobian matrices of the matrix $F\left(\mathbf{x}\right)$ and the matrix $V\left(\mathbf{x}\right)$ at the disease-free equilibrium point, $E_{dfe}$. Thus,

$$J_F\left(E_{dfe}\right)=\left[\begin{array}{cccccc}0& 0& -A_{0}p\left(q-1\right)& 0& -A_0{\eta }_{a}p\left(q-1\right)& 0\\ 0& 0& A_{0}pq& 0& A_0{\eta }_{a}pq& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$$

and

$$J_V\left(E_{dfe}\right)=\left[\begin{array}{cccccc}{A}_3& 0& 0& 0& 0& 0\\ -{\alpha }_u& A_4& 0& 0& 0& 0\\ -{\delta }_u{\sigma }_u& 0& A_5& 0& 0& 0\\ 0& -\gamma {\sigma }_q& -{\varphi }_s& A_6& -{\sigma }_a& 0\\ -A_7& {\sigma }_q\left(\gamma -1\right)& 0& 0& A_8& 0\\ 0& 0& 0& -v_i& 0& A_9\end{array}\right].$$

Therefore, the basic reproduction number, ${\mathcal{R}}_0$, of the model (1) is denoted by:

$$\begin{array}{cc}\hfill {\mathcal{R}}_0& =\frac{p{A}_0}{A_3{A}_4{A}_5{A}_8}((1-\gamma )A_5{\eta }_a{\sigma }_q(A_{3}q-{\alpha }_{u}q+{\alpha }_u)\hfill \\ & +A_4(1-q)(A_5{A}_7{\eta }_a+A_8{\delta }_u{\sigma }_u)).\hfill \end{array}$$

(9)

Theorem 3.

The endemic equilibrium is locally stable in the feasible region (1) if ${\mathcal{R}}_0>1$ [35].

The theorem above establishes a connection between the basic reproduction number and the model’s convergence.

Next, the methodology for controlling the disease, model predictive control (MPC), is discussed, as well as linearized and discretized methods.

3. Model Predictive Control (MPC)

In this section, we look into model predictive control (MPC), which uses the optimum control theory method to obtain the optimal control variables for minimizing the cost function and predicting future output over a time horizon. In addition, the state-space model is generated in the application, and linearization and discretization methods will be used to provide it.

3.1. Linearization and Discretization

The linearization method associated with a Taylor expansion and the exact discretization technique are discussed. The first method is used to convert the model to a linear model, while the second method is used to convert it from a continuous-time model to a discrete-time model.

First, the mathematical model (1) given in continuous time by a system of nonlinear differential equations can be rewritten as

$$\begin{array}{c}\hfill \dot{\mathbf{x}}\left(t\right)=\mathbf{f}(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right)),\end{array}$$

(10)

where $\mathbf{x}$ is the vector of system state, $\mathbf{x}\left(t\right)={\left[S_u\left(t\right)S_{sq}\left(t\right)E_u\left(t\right)E_q\left(t\right)I_s\left(t\right)I_i\left(t\right)I_a\left(t\right)I_{icu}\left(t\right)R\left(t\right)\right]}^T$, $\mathbf{u}\left(t\right)$ is the vector of control input, $\mathbf{u}\left(t\right)={\left[u_{so}\left(t\right)u_M\left(t\right)u_v\left(t\right){\sigma }_a\left(t\right){\alpha }_u\left(t\right)q\left(t\right)\right]}^T$, $\dot{\mathbf{x}}\left(t\right)$ is the vector of the derivative of the system state at time t and $\mathbf{f}$ is a function of the state and input variables of the system.

Let $({\mathbf{x}}_s,{\mathbf{u}}_s)$ be an equilibrium point, i.e., $\mathbf{f}({\mathbf{x}}_s,{\mathbf{u}}_s)=0$, the equilibrium points of the model (1), ${\mathbf{x}}_s$ = $E_{ee}$, the endemic equilibrium point. Furthermore, the initial steady state of the control vector is set as ${\mathbf{u}}_s={\left[000{\sigma }_a{\alpha }_{u}q\right]}^T$. Applying Taylor expansion around the equilibrium point [36],

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}\left(t\right)& =\mathbf{f}\left(\mathbf{x}\right(t),\mathbf{u}(t\left)\right)\hfill \\ \hfill & \simeq \mathbf{f}({\mathbf{x}}_s,{\mathbf{u}}_s)+{\left.\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right|}_{({\mathbf{x}}_s,{\mathbf{u}}_s)}(\mathbf{x}\left(t\right)-{\mathbf{x}}_s)+{\left.\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\right|}_{({\mathbf{x}}_s,{\mathbf{u}}_s)}(\mathbf{u}\left(t\right)-{\mathbf{u}}_s)\hfill \\ \hfill & =A_c{\delta }_{\mathbf{x}}\left(t\right)+B_c{\delta }_{\mathbf{u}}\left(t\right),\hfill \\ \hfill \dot{\mathbf{x}}\left(t\right)& =\frac{d\mathbf{x}\left(t\right)}{dt}=\frac{d(\mathbf{x}\left(t\right)-{\mathbf{x}}_s)}{dt}=\frac{d{\delta }_{\mathbf{x}}\left(t\right)}{dt}=\dot{{\delta }_{\mathbf{x}}}\left(t\right).\hfill \end{array}$$

Therefore, the linear model is

$$\begin{array}{c}\hfill \dot{{\delta }_{\mathbf{x}}}\left(t\right)=A_c{\delta }_{\mathbf{x}}\left(t\right)+B_c{\delta }_{\mathbf{u}}\left(t\right),\end{array}$$

(11)

where $A_c={\left.\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right|}_{({\mathbf{x}}_s,{\mathbf{u}}_s)}$, $B_c={\left.\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\right|}_{({\mathbf{x}}_s,{\mathbf{u}}_s)}$, ${\delta }_{\mathbf{x}}\left(t\right)=\mathbf{x}\left(t\right)-{\mathbf{x}}_s$ and ${\delta }_{\mathbf{u}}\left(t\right)=\mathbf{u}\left(t\right)-{\mathbf{u}}_s$.

Next, the linear model in Equation (11) is discretized by the exact discretization method, which is defined by

$$\begin{array}{c}\hfill {\delta }_{\mathbf{x}}(k+1)=A_d{\delta }_{\mathbf{x}}\left(k\right)+B_d{\delta }_{\mathbf{u}}\left(k\right),\end{array}$$

(12)

where $T_s$ is the sampling period of time, k is discrete-time, $A_d=e^{A_c{T}_s}$ and $B_d=\left(\underset{0}{\overset{T_s}{\int }}{e}^{A_c\tau }d\tau \right)B_c=A_c^{-1}(e^{A_c{T}_s}-I)B_c$ [32]. Finally, the linearized discrete time state-space model of Equation (12) can be rewritten as

$$\begin{array}{cc}\hfill \mathbf{x}(k+1)& =A_d\mathbf{x}\left(k\right)+B_d\mathbf{u}\left(k\right),\hfill \\ \hfill \mathbf{y}\left(k\right)& =C\mathbf{x}\left(k\right),\hfill \end{array}$$

(13)

where $\mathbf{x}$ is the state vector, $\mathbf{u}$ is the input vector and $\mathbf{y}$ is the output vector, with C the matrix indicating the measure state.

3.2. Cost Function of the Model Predictive Control, $J\left(k\right)$

The cost function of model predictive control for constrained problems is defined by

$$\begin{array}{cc}\hfill J\left(k\right)={[Y\left(k\right)-X_{sp}]}^T{W}_y& [Y\left(k\right)-X_{sp}]+▵U^T\left(k\right)W_u▵U\left(k\right),\hfill \\ \hfill \mathrm{s}.\mathrm{t}.\mathbf{x}(k+1)& =A_d\mathbf{x}\left(k\right)+B_d\mathbf{u}\left(k\right),\hfill \\ \hfill \mathbf{y}\left(k\right)& =C\mathbf{x}\left(k\right),\hfill \\ \hfill {\mathbf{u}}_{min}& \le \mathbf{u}\left(k\right)\le {\mathbf{u}}_{max},\hfill \end{array}$$

(14)

where $▵U\left(k\right)=U\left(k\right)-U(k-1)$, $Q_y=diag\left(w_{y_i}\right),Q_u=diag\left(w_{u_j}\right),W_{\mathbf{y}}=I_{N_p}\otimes Q_y,W_{\mathbf{u}}=I_{N_p}\otimes Q_u$ and $X_{sp}={\mathbf{1}}_{N_P}\otimes {\mathbf{x}}_{sp}^T$. The parameters $n_y$ and $n_u$ are the number of output and control variables, $w_{y_i}$, $(i=1,2,...,n_y)$ is the weight of output variables and $w_{u_i}$, $(i=1,2,...,n_u)$ is the weight of control variables. $W_y$ is the weight matrix of the vector, which is the error between the output and the set-point, $W_u$ is the weight matrix of control vector, ${\mathbf{x}}_{sp}$ is target vector or set-point, ${\mathbf{u}}_{min}$ is the minimum of control vector, ${\mathbf{u}}_{max}$ is the maximum of control vector and ⊗ is Kronecker product [32].

The cost function, $J\left(k\right)$, can be transformed into a quadratic programming form, as defined by the equation below:

$$\begin{array}{cc}\hfill \underset{U\left(k\right)}{min}J\left(k\right)& =\frac{1}{2}{U}^T\left(k\right)H_{qp}U\left(k\right)+U^T\left(k\right)G_{qp},\hfill \\ \hfill & H_{qp}=H_y+H_{▵\mathbf{u}},G_{qp}=G_y+G_{▵\mathbf{u}},\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.A_{eq}U\left(k\right)\ge B_{eq},\hfill \end{array}$$

(15)

where $H_y=2{\varphi }^T{W}_x\varphi ,H_{▵\mathbf{u}}=2H_{Q_u},G_y=2({\varphi }^T{W}_{x}F\mathbf{x}(k-1)-{\varphi }^T{W}_x{X}_{sp}),G_{▵\mathbf{u}}=M_u\mathbf{u}(k-1),$$U_{min}={\mathbf{1}}_{N_P}\otimes {\mathbf{u}}_{min}$ and $U_{max}={\mathbf{1}}_{N_P}\otimes {\mathbf{u}}_{max}$. The matrix $A_{eq}$ and the matrix $B_{eq}$ are defined by

$$\begin{array}{c}\hfill A_{eq}=\left[\begin{array}{ccccc}-{\mathbf{I}}_m& \mathbf{O}& \cdots & \mathbf{O}& \mathbf{O}\\ {\mathbf{I}}_m& \mathbf{O}& \cdots & \mathbf{O}& \mathbf{O}\\ \mathbf{O}& -{\mathbf{I}}_m& \cdots & \mathbf{O}& \mathbf{O}\\ \mathbf{O}& {\mathbf{I}}_m& \cdots & \mathbf{O}& \mathbf{O}\\ ⋮& & & & ⋮\\ \mathbf{O}& \mathbf{O}& \cdots & -{\mathbf{I}}_m& \mathbf{O}\\ \mathbf{O}& \mathbf{O}& \cdots & {\mathbf{I}}_m& \mathbf{O}\\ \mathbf{O}& \mathbf{O}& \cdots & \mathbf{O}& -{\mathbf{I}}_m\\ \mathbf{O}& \mathbf{O}& \cdots & \mathbf{O}& {\mathbf{I}}_m\end{array}\right]\mathrm{and}{B}_{eq}=\left[\begin{array}{c}-{\mathbf{u}}_{max}\\ {\mathbf{u}}_{min}\\ -{\mathbf{u}}_{max}\\ {\mathbf{u}}_{min}\\ ⋮\\ -{\mathbf{u}}_{max}\\ {\mathbf{u}}_{min}\\ -{\mathbf{u}}_{max}\\ {\mathbf{u}}_{min}\end{array}\right],\end{array}$$

where ${\mathbf{I}}_m$ is an identity matrix dimension m and $\mathbf{O}$ is a zero matrix dimension m when m is dimension of control vector $\mathbf{u}$.

Next, the COVID-19 epidemic in Thailand is numerically simulated. In addition, the model is subjected to the application of control theory, namely model predictive control.

4. Results

In this section, the baseline parameter values reported in

Table 1. The values of parameters in the model (1).

Parameter	Value	Source	Parameter	Value	Source
${\beta }_{so}$	0.9519	Fitted [11]	$r_s$	1/10	[4]
p	0.8891	Fitted [11]	$r_i$	1/8	[4]
${\sigma }_a$	0.2533	Fitted [11]	$r_a$	0.13978	[4]
${\alpha }_u$	0.0326	Fitted [11]	$r_{icu}$	1/10	[4]
q	0.0397	Fitted [11]	$m_s$	0.00011	Fitted [11]
${\epsilon }_v$	0.75	[37]	$m_i$	0.00011	Fitted [11]
$\mu$	0.000023	[38]	$m_a$	0.5$m_s$	[4]
N	66,190,000	[38]	$m_{icu}$	1.5$m_s$	[4]
$\mathsf{\Lambda }$	$\mu N$		${\varphi }_s$	1/5	[4]
${\psi }_q$	1/14	[4]	$\gamma$	0.7	[4]
${\sigma }_u$	1/5.2	[17]	$v_i$	0.025	[39]
${\sigma }_q$	1/5.2	[17]	${\eta }_a$	0.5	[4]
${\delta }_u$	0.2	[40]	${\mathcal{R}}_0$	1.23

and the initial conditions ($IC$), $S_u\left(0\right)=$ 66,189,981, $S_{sq}\left(0\right)=0$, $E_u\left(0\right)=15$, $E_q\left(0\right)=0$, $I_s\left(0\right)=1$, $I_i\left(0\right)=0$, $I_a\left(0\right)=3$, $I_{icu}\left(0\right)=0$ and $R\left(0\right)=0$, are utilized to fit the model parameters to real data from the COVID-19 pandemic in Thailand. In order to forecast future disease spread behavior, the model parameters are fitted using the sum square error method and the parameters, ${\beta }_{so},p,{\sigma }_a,{\alpha }_u$, and q are computed. Moreover, the model predictive control technique with quadratic programming is used to determine the schedule of COVID-19 disease prevention interventions. The weight matrix of social isolation and vaccination policies are assumed to be 0.01, and the prediction horizon, $N_p$, is 35.

4.1. Mathematical Modeling of COVID-19 Pandemic in Thailand

The COVID-19 pandemic in Thailand is shown in Figure 2. The solid line displays the model’s curve fitting, whereas the dotted line depicts real data during the COVID-19 pandemic in 290 days, from 1 November 2020 to 19 August 2021. The data of newly infected cases in Figure 2a and death cases in Figure 2b are reported by Thailand’s Ministry of Public Health’s Department of Disease Control [11]. In addition, the sum square error method is used to evaluate several parameters in the model (1) relevant to the scenario of disease transmission in Thailand. Thus, parameter values, ${\beta }_{so}=0.9519$, $p=0.8891$, ${\sigma }_a=0.2533$, ${\alpha }_u=0.0326$ and $q=0.0397$, are obtained. Moreover, the basic reproduction number, ${\mathcal{R}}_0$, is calculated as 1.23.

The model parameters were obtained. So, the pandemic trajectory was predicted for 28 weeks later. Figure 3 presents the prediction of infected populations with and without measures, contact tracing of asymptomatic infected and non-quarantined exposed people, as well as the quarantine of non-quarantined susceptible people in the black lines and the red lines. The rise in infected populations is depicted in Figure 3a, with the biggest peak in the number of isolated infected populations, which can affect hospital bed shortages. Furthermore, as time passes, the graph will continue to fluctuate and converge to the endemic equilibrium point. Hence, the results satisfy Theorem 3 for ${\mathcal{R}}_0>1$. In addition, the effectiveness of preventive measures is shown in Figure 3b. The number of infected populations in the graph without control is higher than in the graph with control. Hence, a lot of people will be infected and die because of the pandemic.

4.2. ICU Capacity Restriction and the Effectiveness of Social Isolation Strategies

The results of some components from Figure 3 may be contradictory with the true situation, particularly the number of patients in critical care units (ICUs). The number of ICU admissions is influenced by ICU bed capacity due to ICU bed limits. This means that some patients will have a significant risk of death if they are not treated in an ICU room. In addition, if there are fewer options for increasing the number of ICU beds, the number of serious patient cases must be reduced. Assume that there are roughly 1920 ICU beds available for COVID-19 patients, and set it as a set-point [41]. The model predictive control is then used to create a plan for social isolation regulations in order to keep the number of intensive care unit patients under control.

The 28 weeks of outbreak prediction of infected populations is examined in Figure 4. The solid line in Figure 4a depicts the curve fitting of ICU cases in Thailand. The dashed line depicts the number of ICU admissions, which is influenced by the red line’s capacity for ICU beds. The number of patients is lowered until it is near to the number of ICU beds available, at which point social isolation measures are used to manage the oscillation and convergence to the set-point. As demonstrated in Figure 4b, this result in Figure 4a has an impact on the number of other infected populations.

Figure 5 shows the strategy for preventive measures without vaccination for each week. So, the vaccination rate is set at zero in Figure 5c. The ranges of social distancing and the use of face masks are set between zero and one, while the ranges of ${\sigma }_a$, ${\alpha }_u$, and q are set in zero and their estimated values are shown in Table 1. Moreover, the plan is shown under the relaxation of contact tracing and quarantine. The results pay attention to social distancing and the use of face masks for controlling disease transmission.

Despite the fact that a reduction in the number of patients in ICU classes and all infected people increases the chances of survival, the risk of mortality and the number of infected people remain high. As a result, the effect of managing the epidemic by reducing the number of affected people is investigated.

4.3. The Effectiveness of the Social Isolation Policy in Reducing COVID-19 Pandemic Transmission

Figure 6 illustrates the results of the numbers of infected, susceptible, and recovered populations in Figure 6a,b,c, respectively. The solid lines represent the behavior of a growing population when it is out of control. The dashed lines, on the other hand, depict a population reduction due to control measures. In addition, the effectiveness of social isolation strategies in disease control has led to a reduction in the numbers of infected populations, symptomatic infected, isolated infected, and asymptomatic infected populations, as well as of people in the intensive care unit. The trends of all graphs decrease and approach zero. Moreover, the results of susceptible and recovered populations show a decreasing risk of infection among people.

The measuring plans for each of the two weeks of social isolation strategies without vaccination are shown in Figure 7. The proportions of all social isolation policies are bound by one. The figure suggests the plan for the guidelines aimed at strictly controlling the disease. The plan implies that the disease spread of the COVID-19 pandemic will be controlled if people continue to practice social isolation measures for the rest of their lives. Hence, the policies are good to use for a short time but are unsuitable for dealing with the epidemic over a lengthy period of time.

4.4. The Effectiveness of Social Isolation and Vaccination to Control the COVID-19 Pandemic

We now expect that all Thai citizens will be given immunizations to prevent the disease from spreading further. This study employed the AstraZeneca vaccine, which has a 75 percent efficacy rate and is most likely to be obtained by citizens. The number of all components in the model, except the recovery populations, is kept close to zero in Figure 8. These findings are depicted in Figure 8a,b. The numbers of populations without and under control are represented by the solid and the dashed lines, respectively. The number of sick people will drop and converge to zero as the outbreak is controlled through social isolation and vaccination policies, while the recovery populations will rise. In addition, the number of susceptible populations will decrease. The risk of infection is minimized, and people who are vaccinated will be protected from the COVID-19 pandemic.

The strategy of preventative measures, including vaccination, is shown in Figure 9. Figure 9c suggests and depicts the proportion of people who receive vaccines. As demonstrated in Figure 9a,b, the average proportion of social distancing and the use of face masks are high at first and then decrease. In Figure 9d–f, quarantine and contact tracing are also employed to prevent disease transmission from infected individuals. Moreover, a sufficient quantity of vaccinations and efficient management can improve the situation and quickly bring the pandemic under control.

5. Discussion and Conclusions

This study developed a mathematical model to describe the COVID-19 pandemic’s disease spread in Thailand. The model is built using a system of nonlinear differential equations, and it is then evaluated using model parameter estimation and real data from newly infected and death cases of the COVID-19 pandemic in Thailand. In order to determine some parameters in the model for examining disease spread in Thailand, forecasting future disease is spread to 28 weeks, and analyzing the model, the sum square error method is used to compute the model parameters. The model of the COVID-19 pandemic in Thailand is presented in (1) with the parameters in Table 1. The result of predicting the COVID-19 pandemic in Thailand shows the trend of infected populations in the model and the effectiveness of contact tracing and quarantine. After that, we propose an application of mathematical methods to plan the policy schedule for preventing disease transmission over the time horizon. To control the outbreak, model predictive control is employed to construct a scheme for isolation and vaccination measures for 28 weeks. However, reducing infections completely for all populations in a country is difficult, as it usually depends mostly on the situation in that country, including in Thailand. So, this study investigates the outbreak situation, which includes the limitation of intensive care unit capacity, social isolation without vaccination, and vaccination for a rise in herd immunity in populations.

The first scenario looked at the intensive care unit’s capacity. The number of ICU beds determines the number of patients in intensive care units, and the number of other infected people is decreasing. As a result, each week, policies were able to be relaxed. This suggests that the disease spreads based on the assumption that all patients would be admitted to an ICU service room in order to live. Following that, social isolation strategies were employed to reduce the number of sick people, which caused the number of infected people to drop to near zero. The policy plan is computed every two weeks, and it is suitable to temporarily improve the infection problem. Furthermore, the effectiveness of vaccination is evaluated in the last scenario of this work, and AstraZeneca is used as an example of a vaccine in this study. COVID-19 infection is reduced as a result of social isolation and immunization. The population that receives the vaccination will then be immune to the virus, lowering the risk of illness. The research calculated the social isolation strategy and immunization measures for each of the two weeks. As a result, if this plan can be implemented, disease transmission will be controlled within 28 weeks. As a result, the social isolation police are appropriately utilized to prevent the virus from spreading in a short period of time, but they cannot actually prevent it from spreading because it leads to subsequent difficulties, including economic problems. Vaccination is so critical in the fight against the COVID-19 pandemic.

In conclusion, this study uses mathematical modeling to investigate the COVID-19 pandemic in Thailand. Furthermore, developing preventive strategies over a long time horizon in order to minimize infection in the country depends on the circumstances that may arise. This approach impacts people who can systematically and effectively protect themselves from the virus, which is appropriate for a large group of people to apply in their real life. Furthermore, a proper plan can be implemented to decrease the waste that may result from the implementation of policies.