Adv Syst Sci Appl 2023; 03:66–90 Published online at https://ijassa.ipu.ru. COVID-19 Spread Modeling Incorporating Suggestive Optimal Control Strategies under Uncertainty P.K. Santra1*, D. Pal2, G.S. Mahapatra3, H. Alrabaiah4,5 1Abada Nsup School, Howrah, West Bengal, India 2Chandrahati Dilip Kumar High School, Chandrahati, West Bengal, India 3Department of Mathematics, National Institute of Technology Puducherry, Karaikal, India 4College of Engineering, Al Ain University, Al Ain, United Arab Emirates 5Department of Mathematics, Tafila Technical University, Tafila, Jordan Abstract: In the present paper, we have provided a five-compartmental epidemic model in an interval environment to analyze the spread of COVID-19 infection in India. The proposed model divides the entire population of India into five classes. They are susceptible, exposed, asymptomatic, symptomatic, and recovered classes. Under some suppositions, the crisp model is constructed and converted to an imprecise model by the interval number. We introduced a parametric functional form of an interval number to study the imprecise epidemiological model. The main objective of this study is to develop an epidemiological model in an imprecise environment and to try to understand the dynamics of the epidemic model of COVID-19 infection spread in India. We also presented the COVID-19 model with two controls to effectively control COVID-19 disease in India. Finally, a numerical simulation is carried out considering that the model parameters are imprecise. The numerical results show that our proposed imprecise model is reliable from a practical point of view. Keywords: COVID-19 infection, Basic Reproduction Number, Stability, optimal control, interval number 1. INTRODUCTION The objective of mathematical epidemiology is to find the factors behind the occurrence of a disease. Today, mathematical modeling has played a vital role in understanding epidemiological prototypes of diseases. It also predicts the consequences of the foreword of public health interventions to manage the spread of diseases. In recent times, many researchers and scientists [1–3] have attracted attention to epidemiological systems and have paid great attention to their corresponding research. In an epidemiological system, the transmission of disease from an infected person to a susceptible individual occurs by a mechanism known as the incidence rate, that is, a function describing the above mechanism through different forms. Among them, homogeneous transmission is widely used. This transmission follows the law of mass action. It is also notable that the cost of the treatment strategy to control epidemic diseases is very high. In some cases, proper treatment is only available and inexpensive. Therefore, it is very important to provide the most cost- effective control strategies to prevent the spread of epidemic diseases. Several researchers and scientists [4–6] are trying their best to provide appropriate optimal control strategies for ∗Corresponding author: prasunsantra5@gmail.com COVID-19 SPREAD MODELING UNDER UNCERTAINTY 67 epidemic diseases to protect the human community. Their efforts predict some cost-effective strategies to restrict the extent of epidemic diseases and minimize the cost of the control program. Recently, coronavirus turned out to be a pandemic disease spreading throughout the world and created a panic situation in many countries. COVID-19 infection is a large family of viruses that cause disease like ordinary cold in the primary stage. The latter can create a big problem, such as acute respiratory syndrome. This virus also showed the ability to cause serious health problems among a certain group of individuals, including the elderly population and patients with cardiovascular disease and diabetes. However, the nature of this epidemiology is still changing [7]. On 11 February 2020, this disease was termed coronavirus disease 2019 (COVID-19). The virus behind this disease was declared severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [8]. This virus has been shown to transform from human to human [9]; as a result, the COVID-19 disease spreads day by day around the world [10–14]. Recently, India has become one of the hardest affected country in Asia with COVID- 19 endemic due to its very high population density. The number of COVID positive cases increases progressively, the Indian Government has monitor the succession of this outburst and have taken different public health actions together with social distancing procedures in real time. The figure of positive COVID-19 infection started to increase from 4th March, 2020. As on 8th May, 2020, a total of 59690 confirmed COVID cases together with 17887 recovered and 1986 deaths in India [15]. The Indian government has taken different precaution measures [16–18] to maintain the social distance [19] among the large population of India. To date, there are no specific drugs for COVID-19 infection. Mostly, doctors recommended different treatments via medications to COVID-19 patients depending on their symptoms. As it is not possible to invented anti-COVID therapeutic treatment, it should be restricted via apposite precautionary measures like quarantine mechanisms [21,22], individual safeguard from the infected individual by using social distancing [20], etc. As the COVID-19 virus spread very quickly throughout the world, so various mathematical models depending on the COVID-19 outbreak [23–27] have been performed already. Recntly, Cadoni and Gaeta [28] presented COVID-19 infection in SIR framework. Zhang et al. [29] developed COVID-19 dynamics via fractional derivative. Higazy [30] described COVID-19 virulent disease by using SIDARTHE model. Wu et al. [31] studied a SEIR model to understand the dynamics behind the spread of COVID-19 virus world wise. Read et al. [32] developed a COVID-19 SEIR model based on Poisson-distributed daily time augmentation. Pal et al. [33] explored a COVID-19 based SEQIR model to understand the situation of disease in India. Since the COVID-19 virus has sinister harmful health upshots and pessimistic social impact, a remarkable quantity of awareness and research work are desired on COVID-19 infection. Recently, many researchers [34–36] concentrate their research to develop epidemic models in uncertain / imprecise environments. From their observations, they feel that in reality data cannot be collected precisely due to some unavoidable reasons such as fluctuation of environment, changes in the nature of the virus, ambiguous death rate (like co-morbidity death), uncertainty about the number of undetected infections, human errors, availability of proper information, etc. This paper introduces a five-compartmental COVID-19 infection model in interval environment. We consider parameters of the model that are imprecise in nature and represented by interval numbers. First we formulate the crisp model then we formulate the imprecise model by considering the parameters are interval in nature. We have also introduced two treatment controls in our proposed model under impriseness. The entire paper is organized into several sections and subsections. Section 1 presents the introduction to the paper, which includes a description of the problem and its origin. Section 2 presents the material & Methods which are divided into model variables and parameters, some useful basic mathematical concepts, and the mathematical formulation of the model. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 68 P.K. SANTRA ET AL. The model derivation and preliminaries are explained in Section 2. The basic properties of our proposed model structure are discussed in Section 3. In Section 4, we introduce the concept of the basic reproduction number R0 [37]. Next, we are dealing with two types of equilibrium points in the system, namely disease-free equilibrium (DFE) E0, and endemic equilibrium E1. It is clear that COVID-19 infection is not only a community health trouble [38] but also a great social and monetary shock for developing countries like India. Therefore, it is a very essential concern to control [39–42] the spread of COVID-19 infection in India by adopting an optimal control policy. In Section 5, we have formulated the COVID-19 epidemic model with control treatment. This section provides a procedure for finding an optimal control [43] u(t) that increases the recovery rate and minimizes the cost associated with treatment. The analytical results obtained in the previous sections are numerically verified in Section 6 with the help of realistic values of the model parameters using MATLAB. Lastly, a general conclusion about our proposed model structure is provided in Section 7. 2. MATERIALS & METHODS 2.1. Variables and parameters involved in the Model Following variables and parameters are used to develop the COVID-19 model: Variable Description N(t) : Total populations at time t S(t) : Susceptible people at the time t E(t) : Exposed populations at the time t IA(t) : Infective people in the asymptomatic phase at the time t IS(t) : Infective populations in the symptomatic phase at the time t R(t) : Recovered populations at the time t Parameter Description Λ : Recruitment rate. α : Transmission coefficient due to asymptomatic class (IA). β : Transmission coefficient due to symptomatic class (IS). µ : Natural death rate. ρ : Incubation period. (1− σ) : Fraction of symptomatically infected people (IS). σ : Fraction of asymptotically infected people (IA). (1− k)r : Transfer rate from the asymptomatic phase (IA) to symptomatic phase (IS). kr : Recover rate from the asymptomatic class (IA). m : Recover rate from symptomatic class (IS). d : Death rate from COVID-19 in symptomatic class (IS). Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 69 2.2. Some Useful Definitions Definition 2.1: Interval number: A closed interval number A = [al, ar] in the set of real numbers (R) is defined by A = [al, ar] = {x : al ≤ x ≤ ar, x ∈ R} , where al, and ar are the left and right bound of the interval, respectively. Definition 2.2: Interval Valued Function: The interval valued function [44, 45] of [al, ar] can be expressed as f(p) = a1−p l apr for p ∈ [0, 1]. Arithmetic operations on interval numbers using the concept of interval valued functions are given in Appendix C. 2.3. COVID-19 pandemic model formulation To develop a mathematical model of COVID-19 infection spread in India. The following assumptions and notations are considered throughout the paper. The total human population of India is split up into five mutually exclusive compartments, specifically, susceptible class (S), exposed class (E), infective class in asymptomatic phase (IA), infective class in symptomatic phase (IS) and recovered class (R). Therefore, total population of India N(t) = S(t) + E(t) + IS(t) + IA(t) +R(t). Assumptions: 1. The susceptible population (S) is collected of non-infected persons by COVID-19 still now but can be infected via get in touch with both types of infective (asymptomatic (IA) as well as symptomatic (IS) phase). 2. The exposed population (E) is collected of infected persons by COVID-19 but not infectious. 3. The asymptomatic phase (IA) population is collected of COVID-19 infection persons without any symptoms (but capable of infecting). 4. The infective symptomatic (IS) population is collected of persons with difficulties and diverse symptoms of COVID-19 and need treatment. 5. The recovered (R) population from the disease by treatment or otherwise. 6. The susceptible persons become infected by ample contact with contagious (asymptomatic or symptomatic) persons, and enter into the exposed class. 7. One fraction of asymptomatic infective goes into symptomatic phase and one portion becomes recovered. 8. One fraction of symptomatic contagious goes into recovered class and other portion dies out due to disease related death. 9. All class of population is decreased by natural death. The flow of individuals from one class to another class is presented in Figure 2.1. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 70 P.K. SANTRA ET AL. Fig. 2.1. Flow diagram of individuals of COVID-19 infection 2.4. Crisp model The proposed COVID-19 infection is presented in a crisp environment based on the above stated assumptions through the following set of ordinary differential equations: dS dt = Λ− (αIA + βIS)S − µS dE dt = (αIA + βIS)S − (ρ+ µ)E dIA dt = σρE − (r + µ)IA (2.1) dIS dt = (1− σ)ρE + (1− k)rIA − (m+ µ+ d)IS dR dt = krIA +mIS − µR with their initial conditions : S(0) > 0;E(0) ≥ 0; IA(0) ≥ 0; IS(0) ≥ 0;R(0) ≥ 0. (2.2) 2.5. Imprecise model with interval coefficient Most of the models have been considered in the crips environment, but the data may not be recorded or collected exactly owing to several occasions in reality. Therefore, the model parameters are not to be considered as constants. However, it may be considered as imprecise in nature. To handle with this type of impreciseness, we take the parameters as interval numbers. Therefore, for an imprecise coefficient, we present the epidemic model with an interval counterpart of coefficients, i.e., Λ̂, α̂, β̂, µ̂, ρ̂, r̂, m̂ and d̂. Then the COVID-19 Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 71 infection model (2.1) with interval valued parameters modified to the following form: dS dt = Λ̂− ( α̂IA + β̂IS ) S − µ̂S dE dt = ( α̂IA + β̂IS ) S − (ρ̂+ µ̂)E dIA dt = σρ̂E − (r̂ + µ̂)IA (2.3) dIS dt = (1− σ)ρ̂E + (1− k)r̂IA − (m̂+ µ̂+ d̂)IS dR dt = kr̂IA + m̂IS − µ̂R with their initial conditions of equation (2.2), where Λ̂ ∈ [Λl,Λu], α̂ ∈ [αl, αu], β̂ ∈ [βl, βu], µ̂ ∈ [µl, µu], ρ̂ ∈ [ρl, ρu], r̂ ∈ [rl, ru], m̂ ∈ [ml,mu] and d̂ ∈ [dl, du]. Also Λl > 0, αl > 0, βl > 0, µl > 0, ρl > 0, rl > 0, ml > 0, dl > 0. Using the parametric functional form of the interval number, the pandemic system (2.3) in the parametric form is as follows: dS dt = Λ1−p l Λp u − { α1−p l αp uIA + β1−p l βp uIS } S − µ1−p l µp uS dE dt = { α1−p l αp uIA + β1−p l βp uIS } S − { ρ1−p l ρpu + µ1−p l µp u } E dIA dt = σρ1−p l ρpuE − { r1−p l rpu + µ1−p l µp u } IA (2.4) dIS dt = (1− σ)ρ1−p l ρpuE + (1− k)r1−p l rpuIA − { m1−p l mp u + µ1−p l µp u + d1−p l dpu } IS dR dt = kr1−p l rpuIA +m1−p l mp uIS − µ1−p l µp uR for p ∈ [0,1]. The above system formulation (2.4) can be rewritten as: dS dt = Λ1−p l Λp u − { α1−p l αp uIA + β1−p l βp uIS } S − µ1−p l µp uS dE dt = { α1−p l αp uIA + β1−p l βp uIS } S − AE dIA dt = σρ1−p l ρpuE −BIA (2.5) dIS dt = (1− σ)ρ1−p l ρpuE + (1− k)r1−p l rpuIA − CIS dR dt = kr1−p l rpuIA +m1−p l mp uIS − µ1−p l µp uR Where A = [ ρ1−p l ρpu + µ1−p l µp u ] ; B = [ r1−p l rpu + µ1−p l µu) p ] ; C = [(ml) 1−p(mu) p + (µl) 1−p(µu) p + (dl) 1−p(du) p] and p ∈ [0,1], with initial conditions as stated in equation (2.2). Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 72 P.K. SANTRA ET AL. 3. THEORETICAL STUDY OF THE MODEL 3.1. Positivity of the solutions Theorem 3.1: Each solution of the system (2.5) under conditions (2.2) exists in the interval [0,∞) and initial conditions (2.2) are satisfied for all t ≥ 0. Proof In view of the fact that the right-hand side of the parametric pandemic system (2.5) is com- pletely continuous and locally Lipschitzian on C, the solution (S(t), E(t), IA(t), IS(t), R(t)) of (2.5) with initial conditions (2.2) exists and is unique on [0, κ), where 0 < κ < +∞. First, we show that S(t) > 0, ∀ t ∈ [0, κ). Otherwise, ∃ a t ∈ (0, κ) such that S(t1) = 0, · S(t1) ≤ 0 and S(t) > 0 for all t ∈ [0, t1). Hence, there must have E(t) ≥ 0 for all t ∈ [0, t1). If this statement is not true, then there exists a t2 ∈ (0, t1) such that E(t2) = 0, · E(t2) < 0 and E(t) ≥ 0 on [0, t2). We claim that IA(t) ≥ 0, ∀ t ∈ [0, t2). If this is not true, then ∃ a t3 ∈ (0, t2) such that IA(t3) = 0, · IA(t3) < 0 and IA(t) ≥ 0 on [0, t3). Now, from the third equation of (2.5), we have dIA(t3) dt = σ(ρl) 1−p(ρu) pE(t3)−BIA(t3) = σ(ρl) 1−p(ρu) pE(t3) ≥ 0, which is a contradiction with · IA(t3) < 0. So IA(t) ≥ 0, ∀ t ∈ [0, t2). Next, we assert that IS(t) ≥ 0, ∀ t ∈ [0, t2) . If this is not correct, then ∃ a t4 ∈ (0, t2) such that IS(t4) = 0, · IS(t4) < 0 and IS(t) ≥ 0 on [0, t4). Now, from the fourth equation of (2.5), we have dIS(t4) dt = (1− σ)(ρl) 1−p(ρu) pE(t4) + (1− k)(rl) 1−p(ru) pIA(t4)− CIS(t4) = (1− σ)(ρl) 1−p(ρu) pE(t4) + (1− k)(rl) 1−p(ru) pIA(t4) ≥ 0, which contradict the consideration · IS(t4) < 0, therefore IS(t) ≥ 0, ∀ t ∈ [0, t2). Now, using the second equation of (2.5), we have dE(t2) dt = ( (αl) 1−p(αu) pIA(t2) + (βl) 1−p(βu) pIS(t2) ) S(t2)− AE(t2) = ( (αl) 1−p(αu) pIA(t2) + (βl) 1−p(βu) pIS(t2) ) S(t2) ≥ 0 which is a contradiction with · E(t2) < 0, therefore E(t) ≥ 0, ∀ t ∈ [0, t1). Hence, E(t), IA(t), IS(t) ≥ 0, ∀ t ∈ [0, t1). Next, we claim that R(t) ≥ 0, ∀ t ∈ [0, t1). If this is not true, then ∃ a t5 ∈ (0, t1) such that R(t5) = 0, · R(t5) < 0 and R(t) ≥ 0 on [0, t5) . From the last equation of (2.5), we have dR(t5) dt = k(rl) 1−p(ru) pIA(t5) + (ml) 1−p(mu) pIS(t5)− (µl) 1−p(µu) pR(t5) = k(rl) 1−p(ru) pIA(t5) + (ml) 1−p(mu) pIS(t5) ≥ 0, which contradict the consideration · R(t5) < 0, therefore R(t) ≥ 0, ∀ t ∈ [0, t1). Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 73 Again, from the first equation of (2.5) follows that dS(t1) dt = Λ1−p l Λp u − ( α1−p l αp uIA(t1) + β1−p l βp uIS(t1) ) S(t1)− µ1−p l µp uS(t1) = Λ1−p l Λp u > 0, This shows that S(t) > 0, ∀ t ∈ [0, κ). Based on the steps discussed previously, it observes that E(t) ≥ 0, IA(t) ≥ 0, IS(t) ≥ 0, R(t) ≥ 0 for all t ∈ [0, κ), where 0 < κ < +∞. This completes the proof. 3.2. Invariant region Theorem 3.2: The feasible region µ defined by µ = { (S(t), E(t), IA(t), IS(t), R(t)) ∈ R6 + : 0 < N ≤ Λ1−p l Λp u ϑ } where ϑ = min {(µl) 1−p(µu) p, (µl) 1−p(µu) p + (dl) 1−p(du) p}, with intial conditions (2.2) is positively invariant. Proof By adding the equations of the system (2.5) we obtain dN dt = (Λl) 1−p(Λu) p − (µl) 1−p(µu) p(S + E + IA + IS +R)− (dl) 1−p(du) pIS Therefore, dN dt + ϑN = Λ1−p l Λp u − (µ1−p l µp u − ϑ)(S + E + IA +R)− (µ1−p l µp u + d1−p l dpu − ϑ)IS ≤ (Λl) 1−p(Λu) p (3.6) where ϑ = min {(µl) 1−p(µu) p, (µl) 1−p(µu) p + (dl) 1−p(du) p}. The solution N(t) of the differential equation (3.6) has the following property, 0 < N(t) ≤ N(0) exp(−ϑt) + (Λl) 1−p(Λu) p ϑ (1− exp(−ϑt)), where N(0) represents the sum of the initial values of the variables. As t → ∞, we have 0 < N(t) ≤ (Λl) 1−p(Λu)p ϑ . Moreover, if N(0) ≤ (Λl) 1−p(Λu)p ϑ then also N(t) ≤ (Λl) 1−p(Λu)p ϑ for all t. This means that (Λl) 1−p(Λu)p ϑ is the upper bound of N . On the other hand, if N(0) > (Λl) 1−p(Λu)p ϑ , then N(t) will decrease to (Λl) 1−p(Λu)p ϑ . This means that if N(0) > (Λl) 1−p(Λu)p ϑ , then the solution (S(t), E(t), IA(t), Is(t), R(t)) enters µ or approaches it asymptotically. Hence, it is positively invariant under the flow induced by the systems (2.5) with initial conditions (2.2). Thus in µ the model (2.5) with initial conditions (2.2) is well-posed both mathematically and epidemiologically. Therefore, the study of the dynamics of the proposed COVID-19 system in µ is mathematically sufficient. 3.3. Existence of Equilibria and Basic reproduction number We will study the existence and stability behavior of the system (2.5) at different equilibrium points. The equilibrium points of the proposed system are as follows: Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 74 P.K. SANTRA ET AL. (i) Disease-free equilibrium (DFE): E0( (Λl) 1−p(Λu)p (µl)1−p(µu)p , 0, 0, 0, 0), (ii) Endemic equilibrium: E1(S ∗, E∗, I∗A, I ∗ S, R ∗). The Basic Reproduction Number (BRN): The basic reproduction number [46] is defined as ”the number of new infective individuals produced by a single infective individual during his or her effective infectious period when introduced into susceptible populations”. Now, the BRN of COVID-19 system (2.5) will be derived using the next generation matrix method [47]. Let z = (E(t), IA(t), IS(t), R(t), S(t))T , the proposed pandemic system (2.5) can be written in the following form: dz dt = 𭟋(z)− υ(z) where 𭟋(z) =  ((αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS)S 0 0 0 0  and υ(z) =  AE −σ(ρl) 1−p(ρu) pE +BIA −(1− σ)(ρl) 1−p(ρu) pE − (1− k)(rl) 1−p(ru) pIA + CIS −k(rl) 1−p(ru) pIA − (ml) 1−p(mu) pIS + (µl) 1−p(µu) pR −(Λl) 1−p(Λu) p + ((αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS + (µl) 1−p(µu) p)S  The Jacobian matrices of 𭟋(z) and υ(z) at E0 are as follows: D𭟋(E0) = [ F3×3 0 0 0 0 0 0 0 0 ] and Dυ(E0) =  V3×3 0 0 0 −kr −(ml) 1−p(mu) p (µl) 1−p(µu) p 0 0 (Λl) 1−p(Λu)pα (µl)1−p(µu)p (Λl) 1−p(Λu)pβ (µl)1−p(µu)p 0 (µl) 1−p(µu) p  where F = 0 (Λl) 1−p(Λu)pα (µl)1−p(µu)p (Λl) 1−p(Λu)pβ (µl)1−p(µu)p 0 0 0 0 0 0  and V = [ A 0 0 −σ(ρl) 1−p(ρu) p B 0 −(1− σ)(ρl) 1−p(ρu) p −(1− k)r C ] Now, FV −1 is the next generation matrix of (2.5), consequently the spectral radius of the matrix FV −1 is denoted and defined by [47]: r ( FV −1 ) = Λ1−p l Λp u µ1−p l µp u [α1−p l αp uσρ 1−p l ρpuC + β1−p l βp uρ 1−p l ρpu{σ(1− k)r1−p l rpu +B(1− σ)}] ABC Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 75 According to Ref. [47], the BRN of system (2.5) is R0 = r ( FV −1 ) = Λ1−p l Λp u µ1−p l µp u [α1−p l αp uσρ 1−p l ρpuC + β1−p l βp uρ 1−p l ρpu{σ(1− k)r1−p l rpu +B(1− σ)}] ABC Notice that (Λl) 1−p(Λu)p (µl)1−p(µu)p is the number of susceptibles at the DFE. Existence of Endemic Equilibrium E1(S ∗, E∗, I∗A, I ∗ S, R ∗): Here, we analyze the existence of nontrivial endemic equilibrium E1(S ∗, E∗, I∗A, I ∗ S, R ∗) of the system (2.5). To find the endemic equilibrium of pandemic model (2.5), we consider the following: S(t) > 0, E(t) > 0, IA(t) > 0, IS(t) > 0, R(t) > 0 and dS dt = 0, dE dt = 0, dIA dt = 0, dIS dt = 0, dR dt = 0 (3.7) The solution of the equations of system (3.7) is: S∗ = ABC (αl)1−p(αu)pσ(ρl)1−p(ρu)pC+(βl)1−p(βu)p(ρl)1−p(ρu)p{σ(1−k)(rl)1−p(ru)p+B(1−σ)} = (Λl) 1−p(Λu)p (µl)1−p(µu)p 1 R0 , E∗ = BI∗A σ(ρl)1−p(ρu)p , I∗A = (Λl) 1−p(Λu)pσ(ρl) 1−p(ρu)p(1− 1 R0 ) AB , I∗S = {σ(1−k)(rl) 1−p(ru)p+B(1−σ)}I∗A σC , R∗ = [k(rl) 1−p(ru)pσC+(ml) 1−p(mu)p{σ(1−k)(rl) 1−p(ru)p+B(1−σ)}]I∗A σ(µl)1−p(µu)pC , Hence, I∗A has positive solution iff R0 > 1. Summarizing the above discussions, we arrive at the following result. Theorem 3.3: The system (2.5) has a DFE E0( Λ1−p l Λp u µ1−p l µp u , 0, 0, 0, 0) for all parameter values. If R0 > 1, the system (2.5) also asserts a unique endemic equilibrium E1(S ∗, E∗, I∗A, I ∗ S, R ∗). 3.4. Stability Behavior Here, we examine the condition of stability of the pandemic system (2.5) in two different equilibrium points E0 and E1. For the stability of DFE E0( (Λl) 1−p(Λu)p (µl)1−p(µu)p , 0, 0, 0, 0) we consider the theorems given below Theorem 3.4: The DFE E0 of the system (2.5) is locally asymptotically stable if R0 < 1. Proof Detailed proof of this theorem is given in Appendix A. Theorem 3.5: The DFE E0( (Λl) 1−p(Λu)p (µl)1−p(µu)p , 0, 0, 0, 0) is globally asymptotically stable for the system (2.5) if R0 < 1 and became unstable if R0 > 1. Proof The COVID-19 system (2.5) rewrite as dX dt = F (X, V ) dV dt = G(X, V ), such that G(X, 0) = 0 Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 76 P.K. SANTRA ET AL. where X = (S,R) ∈ R2 denotes compartments of the number of uninfected individuals, V = (E(t), IA(t), IS(t)) ∈ R4 is the number of infected individuals compartments, and E0( (Λl) 1−p(Λu)p (µl)1−p(µu)p , 0, 0, 0, 0) is the DFE of the system (2.5). The global stability of the DFE is guaranteed if satisfying the following two conditions: (i) For dX dt = F (X, 0), X∗ is globally asymptotically stable, (ii) G(X, V ) = BV − Ĝ(X, V ), Ĝ(X, V ) ≥ 0 for (X, V ) ∈ Ω where B = DVG(X∗, 0) is a Metzler matrix and Ω is the positively invariant set of the system(2.5). Following Castillo-Chavez et al. [48], we investigate the aforementioned conditions for pandemic system (2.5), F (X, 0) = [ (Λl) 1−p(Λu) p − (µl) 1−p(µu) pS 0 ] B =  −A (Λl) 1−p(Λu)pα (µl)1−p(µu)p (Λl) 1−p(Λu)pβ (µl)1−p(µu)p σ(ρl) 1−p(ρu) p −B 0 (1− σ)(ρl) 1−p(ρu) p (1− k)(rl) 1−p(ru) p −C  and Ĝ(X, V ) = [ 0 0 0 ] As the off-diagonal elements of B are non negative, and Ĝ(X, V ) ≥ 0, hence the disease-free equilibrium E0( (Λl) 1−p(Λu)p (µl)1−p(µu)p , 0, 0, 0, 0) is globally asymptotically stable if R0 < 1 and became unstable if R0 > 1. Theorem 3.6: The endemic equilibrium point E1(S ∗, E∗, I∗A, I ∗ S, R ∗) of the system (2.5) is locally asymptotically stable if R0 > 1, B1B2 −B3 > 0 and B1B2B3 −B2 1B4 −B2 3 > 0 Proof Detailed proof of this theorem is given in Appendix B. 4. OPTIMAL CONTROL FOR COVID-19 MODEL The most important motive of studying infectious diseases is to control the transmission of the disease and finally stop the infection from the population. transmission of the disease and finally stop the infection from the population. In the present scenario, many developed countries are working hard to discover the medicine or vaccine of COVID-19. As vaccine of COVID-19 has not yet been discovered, we have to control this pandemic by suitable controlling strategies like the use of proper drugs, etc. This section presents an optimal control problem relative to the COVID-19 epidemic model (2.5) with two control variables u1(t) and u2(t), as it is essential to construct an optimal control problem in such a way that the total amount of drugs can be minimized. The first control u1(t) represents the recovery rate of the infectious asymptomatic individuals with treatment. The second control u2(t) represents the recovery rate of the infectious symptomatic individuals with treatment. Therefore, the proposed COVID-19 epidemic model with two controls becomes: Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 77 dS dt = Λ− ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) S − (µl) 1−p(µu) pS dE dt = ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) S − ((ρl) 1−p(ρu) p + (µl) 1−p(µu) p)E dIA dt = σ(ρl) 1−p(ρu) pE − (µ+ (µl) 1−p(µu) p + u1)IA (4.8) dIS dt = (1− σ)(ρl) 1−p(ρu) pE + µIA − (u2 + (µl) 1−p(µu) p + (dl) 1−p(du) p)IS dR dt = u1IA + u2IS − (µl) 1−p(µu) pR satisfying, S(0) = S0, E(0) = E0, IA(0) = IA0 , IS(0) = IS0 , R(0) = R0 (4.9) Here µ stands for the rate of transfer from the asymptomatic phase to the symptomatic phase of the infective individual, and all other system parameters remain unchanged as described earlier. Let µ̂ be the interval counterparts of µ, where µ̂ ∈ [µl, µu]. Parametric form of interval valued parameter µ̂ is (µl) 1−p(µu) p for p ∈ [0,1]. The objective functional [50] is defined as J(u1(t), u2(t)) = tf∫ 0 [G1E + IAE +G3IS + K1 2 u2 1 + K2 2 u2 2]dt (4.10) where G1,G2,G3,K1 2 and K2 2 are positive constants. The square of the control variable reflects the severity of the side effects of the treatment. We have to minimize the objective functional J(u1(t), u2(t)) given in (4.8) so that the COVID infected individuals and the cost of treatment can be minimized. Therefore, we look for an optimal control (u∗ 1, u ∗ 2) such that J(u∗ 1, u ∗ 2) = min {J(u1, u2) : (u1, u2) ∈ U} (4.11) where admissible control set U = {(u1, u2) : uiis measurable, 0 ≤ ui ≤ 1, t ∈ [0, tf ], fori = 1, 2}. 4.1. Existence of an optimal control Theorem 4.1: An optimal control (u∗ 1, u ∗ 2) exists for J(u∗ 1, u ∗ 2) = min {J(u1, u2) : (u1, u2) ∈ U} for the control system (4.8) with initial conditions (4.9). Proof This theorem can be proved by using the results of the existence of an optimal control pair from Ref. [51]. The system of equation (4.8) is bounded from above by a linear system. The boundedness of solutions of system (4.8) for a finite time interval is used to prove the existence of an optimal control. We first check the following properties before using of the theorem from Ref. [51]: (i) Set of state variables and their control is non-empty, (ii) U is convex and closed, (iii) The right-hand side of (4.8) is bounded by a linear function in the state control, (iv) The integrand of the objective functional is convex on U , (v) There exist constants c1, c2 > 0 and q > 1 such that the integrand of (4.10) satisfies G1E +G2IA +G3IS + K1 2 u2 1 + K2 2 u2 2 ≥ c1(u 2 1 + u2 2) q 2 − c2. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 78 P.K. SANTRA ET AL. We use a result by Lukes ( [52] Theorem 9.2.1) to verify the condition (i) for the system (4.8) with bounded coefficients. The control set U is convex and closed by definition, which gives the condition (ii). The right-hand side of the state system (4.8) satisfies condition (iii) as the state solutions are a priori bounded. The integrand in the objective functional, G1E +G2IA +G3IS + K1 2 u2 1 + K2 2 u2 2, is clearly convex on U , which gives condition (iv). Finally, because the state variables are bounded, then there are c1, c2 > 0 and q > 1 satisfying G1E +G2IA +G3IS + K1 2 u2 1 + K2 2 u2 2 ≥ c1(u 2 1 + u2 2) q 2 − c2 Based on the above conditions, there exists an optimal control (u∗ 1, u ∗ 2) such that J(u∗ 1, u ∗ 2) = min {J(u1, u2) : (u1, u2) ∈ U} 4.2. Characterization of the optimal control The Pontryagin’s Maximum Principle [53] is used to derive the necessary conditions for the optimal control pair. The necessary conditions for optimum control are obtained by means of Pontryagin’s Maximum Principle. Theorem 4.2: There exists an optimal control (u∗ 1, u ∗ 2) and corresponding solutions S ∗ , E ∗ , IA ∗ , IS ∗ , R ∗ that minimizes J of (4.10) over U . The explicit optimal controls are connected to the existence of continuous specific functions λi(t), for i = 1, 2, 3, 4, 5, the solutions of the following adjoint system: dλ1 dt = (λ1 − λ2) ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) + λ1(µl) 1−p(µu) p dλ2 dt = −G1 + λ2((ρl) 1−p(ρu) p + (µl) 1−p(µu) p)− (λ3 − λ4)σ(ρl) 1−p(ρu) p − λ4(ρl) 1−p(ρu) p dλ3 dt = −G2 + (λ1 − λ2)α 1−p l αp uS + (λ3 − λ5)u1(t) + (λ3 − λ4)(µl) 1−p(µu) p − λ3(µl) 1−p(µu) p dλ4 dt = −G3 + (λ1 − λ2)(βl) 1−p(βu) pS + (λ4 − λ5)u2(t) + λ4(µl) 1−p(µu) p dλ5 dt = λ5(µl) 1−p(µu) p On the conditions of transversality: λi(tf ) = 0, for i = 1, 2, 3, 4, 5. In addition, the following property holds: u∗ 1 = min{max{0, (λ3 − λ5)IA ∗ K1 }, 1}, and u∗ 2 = min{max{0, (λ4 − λ5)IS ∗ K2 }, 1}. Proof The Hamiltonian is defined as follows: Ĥ = G1E +G2IA +G3IS + K1 2 u2 1 + K2 2 u2 2 + λ1[Λ 1−p l Λp u − (α1−p l αp uIA + β1−p l βp uIS)S(4.12) −µ1−p l µp uS] + λ2[(α 1−p l αp uIA + β1−p l βp uIS)S − (ρ1−p l ρpu + µ1−p l µp u)E] +λ3[σρ 1−p l ρpuE − (µ1−p l µp u + µ1−p l µp u + u1)IA] + λ4[(1− σ)ρ1−p l ρpuE + µ1−p l µp uIA −(u2 + µ1−p l µp u + d1−p l dpu)IS] + λ5[u1IA + u2IS − (µl) 1−p(µu) pR] Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 79 where λi (i = 1, 2, 3, 4, 5, 6) are the adjoint functions to be determined suitably. The form of the equations of adjoint and transversality conditions are the standard consequences from Pontryagin’s Maximum Principle [53]. The adjoint system for the proposed COVID-19 system can be obtained as follows: dλ1 dt = −∂Ĥ ∂S = (λ1 − λ2) ((αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS) + λ1(µl) 1−p(µu) p dλ2 dt = −∂Ĥ ∂E = −G1 + λ2 ( ρ1−p l ρpu + µ1−p l µp u ) − (λ3 − λ4)σ(ρl) 1−p(ρu) p − λ4(ρl) 1−p(ρu) p dλ3 dt = − ∂Ĥ ∂IA = −G2 + (λ1 − λ2)α 1−p l αp uS + (λ3 − λ5)u1(t) + (λ3 − λ4)µ 1−p l µp u + λ3µ 1−p l µp u dλ4 dt = − ∂Ĥ ∂IS = −G3 + (λ1 − λ2)β 1−p l βp uS + (λ4 − λ5)u2(t) + λ4((µl) 1−p(µu) p + (dl) 1−p(du) p) dλ5 dt = −∂Ĥ ∂R = λ5(µl) 1−p(µu) p (4.13) The transversality conditions (or boundary conditions) are λi(tf ) = 0, i = 1, 2, 3, 4, 5 (4.14) Using the optimality condition, we have ∂Ĥ ∂u1 = K1u ∗ 1 − (λ3 − λ5)I ∗ A = 0 at u1 = u∗ 1(t) ⇒ u∗ 1(t) = (λ3−λ5)I ∗ A K1 ∂Ĥ ∂u2 = K2u ∗ 2 − (λ4 − λ5)I ∗ S = 0 at u2 = u∗ 2(t) ⇒ u∗ 2(t) = (λ4−λ5)I ∗ S K2 (4.15) Again, from the bounds for the control u1(t), we get u∗ 1 =  (λ3−λ5)I ∗ A K1 , if 0 ≤ (λ3−λ5)I ∗ A K1 ≤ 1 0, if (λ3−λ5)I ∗ A K1 ≤ 0 1, if (λ3−λ5)I ∗ A K1 ≥ 1 In compact form: u∗ 1 = min{max{0, (λ3 − λ5)I ∗ A K1 }, 1} (4.16) From the bounds for the control u2(t), we get u∗ 2 =  (λ4−λ5)I ∗ S K2 , if 0 ≤ (λ4−λ5)I ∗ S K2 ≤ 1 0, if (λ4−λ5)I ∗ S K2 ≤ 0 1, if (λ4−λ5)I ∗ S K2 ≥ 1 The above can be written in dense form: u∗ 2 = min{max{0, (λ4 − λ5)I ∗ S K2 }, 1} (4.17) Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 80 P.K. SANTRA ET AL. Using (4.16), we obtain the following optimality system: dS dt = (Λl) 1−p(Λu) p − ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) S − (µl) 1−p(µu) pS (4.18) dE dt = ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) S − ( (ρl) 1−p(ρu) p + (µl) 1−p(µu) p ) E dIA dt = σ(ρl) 1−p(ρu) pE − ((µl) 1−p(µu) p + (µl) 1−p(µu) p)IA −min{max{0, (λ3 − λ5)I ∗ A K1 }, 1}IA dIS dt = (1− σ)ρ1−p l ρpuE + µ1−p l µp uIA −min{max{0, (λ4 − λ5)I ∗ S K2 }, 1}IS − (µ1−p l µp u + d1−p l dpu)IS dR dt = min{max{0, (λ3 − λ5)I ∗ A K1 }, 1}IA +min{max{0, (λ4 − λ5)I ∗ S K2 }, 1}IS − (µl) 1−p(µu) pR dλ1 dt = (λ1 − λ2) ( (αl) 1−p(αu) pIA + (βl) 1−p(βu) pIS ) + λ1(µl) 1−p(µu) p dλ2 dt = −G1 + λ2 ( (ρl) 1−p(ρu) p + (µl) 1−p(µu) p ) − (λ3 − λ4)σ(ρl) 1−p(ρu) p − λ4(ρl) 1−p(ρu) p dλ3 dt = −G2 + (λ1 − λ2)α 1−p l αp uS + (λ3 − λ5)min{max{0, (λ3 − λ5)I ∗ A K1 }, 1}+ (λ3 − λ4)µ 1−p l µp u + λ3µ 1−p l µp u dλ4 dt = −G3 + (λ1 − λ2)β 1−p l βp uS + (λ4 − λ5)min{max{0, (λ4 − λ5)I ∗ S K2 }, 1}+ λ4(µ 1−p l µp u + d1−p l dpu) dλ5 dt = λ5(µl) 1−p(µu) p subject to the following conditions: S(0) = S0, E(0) = E0, IA(0) = IA0, IS(0) = IS0, R(0) = R0 and λ1(tf ) = 0, λ2(tf ) = 0, λ3(tf ) = 0, λ4(tf ) = 0, λ5(tf ) = 0 This completes the proof. 5. NUMERICAL SIMULATIONS In this section, we analyze our mathematical model in different environment through some simulation works. The simulation works are done by using MATLAB R2008a software package with the use of the most versatile ordinary differential equation solver ODE45. Simulations of Crips COVID-19 Model: To understand the system dynamics of our proposed crisp model more clearly, we present the following graphical presentation based on Table 5.1 and Table 5.2. Table 5.1. Model parameters for COVID-19 system Parameters Values References Λ 60000 Assumed α 0.0000000001 Fitted β 0.0000000001 Fitted µ 0.00002 Assumed ρ 0.2 Fitted r 0.1 Fitted m 0.07 Fitted d 0.0018 Fitted σ 0.64 Fitted k 0.7 Fitted Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 81 Table 5.2. Preliminary population density for COVID-19 mode S(0) 1200000000 E(0) 200000 IA(0) 200000 IS(0) 220114 R(0) 347979 Using the parameter values given in Table 5.1, the endemic equilibrium point is given by (713600234.09, 228324.60, 291846.24, 350560.38, 2251062473.25) and correspoinding values of R0 is given by 4.1999 > 1. Therefore, all population of our system is all persist, i.e., the population of our model system converges to the corresponding endemic equilibrium point (713600234.09, 228324.60, 291846.24, 350560.38, 2251062473.25). Also, since R0 > 1, indicates the presence of infection in the system. Using the parameter values given in Table 5.1 for different initial conditions given in Table 5.2, the dynamics of the model is presented in Fig. 5.2 and Fig. 5.3. 0 50 100 150 200 250 300 350 400 450 500 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 7 Time A c ti v e i n fe c te d ( I S ) Fig. 5.2. Time series plot of the active infected population (IS) for the parameter values given in Table 1 and the initial conditions given in Table 2. 0 100 200 300 400 500 2 4 6 8 10 12 x 10 8 Time S u s c e p ti b le ( S ) 0 100 200 300 400 500 0 1 2 3 4 x 10 7 Time E x p o s e d ( E ) 0 100 200 300 400 500 0 1 2 3 4 5 x 10 7 Time A s y m p to m a ti c ( I A ) 0 100 200 300 400 500 0 2 4 6 8 10 x 10 8 Time R e c o v e re d ( R ) i ii iviii Fig. 5.3. Time series plot of (i) susceptible population, (ii) exposed population, (iii) infective population asymptomatic phase, (iv) recovered population based on parameter values given in Table 1 and initial conditions given in Table 2. From Fig. 5.2, it is observed that the pick of height is nearly equal to 4.7× 107. Moreover, after reaching the height pick of the infection, the active infective curve is gradually decreased as the time progress. Furthermore, the curve of exposed population (E), the infective population in asymptomatic phase (IA) shows the similar behavior as infective (IS) curve for different with respect to time t (see Fig. 5.3). On the contrary, Susceptible (S) population curve gradually decreases and gradually increases as the time progress. Imprecise Model To verify our analytical findings, we set the values of the parameters as interval numbers given in Table 5.1, and the initial population density is given in Table 5.2, respectively. Based on Table 5.1, we calculate the Disease-free equilibrium, Endemic equilibrium, and R0 for different values of p ∈ [0, 1] which are given in Table 5.3. From Table 5.3, we observe that E0 and E1 exist for all values of p. It is also interesting to observe that the value of R0 is continuously increases with the increasing value of p. It Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 82 P.K. SANTRA ET AL. Table 5.3. Disease-free equilibrium (E0), Endemic equilibrium (E1) and R0 p Disease-free equilibrium Endemic equilibrium R0 0.0 (2947368421.05, 0, 0, 0, 0) (778798079.86, 216835.35, 266284.03, 321556.78, 2138994795.46) 3.7845 0.1 (2957243264.96, 0, 0, 0, 0) (765297598.12, 219172.65, 271318.04, 327289.35, 2161811456.07) 3.8642 0.2 (2967151193.49, 0, 0, 0, 0) (752031614.29, 221489.81, 276390.65, 333055.41, 2184422467.99) 3.9455 0.3 (2977092317.47, 0, 0, 0, 0) (738996051.42, 223787.25, 281502.46, 338855.54, 2206831952.79) 4.0286 0.4 (2987066748.13, 0, 0, 0, 0) (726186903.48, 226065.38, 286654.12, 344690.33, 2229043960.98) 4.1134 0.5 (2997074597.06, 0, 0, 0, 0) (713600234.09, 228324.60, 291846.24, 350560.38, 2251062473.25) 4.1999 0.6 (3007115976.22, 0, 0, 0, 0) (701232175.32, 230565.32, 297079.47, 356466.27, 2272891401.66) 4.2883 0.7 (3017190997.95, 0, 0, 0, 0) (689078926.54, 232787.93, 302354.44, 362408.60, 2294534590.87) 4.3786 0.8 (3027299774.97, 0, 0, 0, 0) (677136753.17, 234992.80, 307671.80, 368387.99, 2315995819.22) 4.4707 0.9 (3037442420.37, 0, 0, 0, 0) (665401985.62, 237180.32, 313032.22, 374405.02, 2337278799.96) 4.5648 1.0 (3047619047.62, 0, 0, 0, 0) (653871018.08, 239350.87, 318436.33, 380460.31, 2358387182.31) 4.6609 is also interesting to see that BRN of our proposed number is gradually increasing with the increasing value of p. Table 5.3 also indicates that BRN R0 is always much greater than 1 for any values of p. Hence, the infection spread very quickly. Therefore, Government should take proper policy to reduce the value of R0 below or near to 1. Again, we see from Table 5.3 that our model system’s population converges to the corresponding endemic equilibrium points E1 for all values of p (by our analytical results stated in Theorem 6). Graphical presentation of BRN R0 with respect to is presented through Fig. 5.4. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 p R 0 Fig. 5.4. Reproduction number R0 with respect to p This figure shows that as the parameter p increases (0 to 1) BRN is also increases i.e., speeding rate of infection are increased as the parameter p increases. Therefore, from Table 5.3 and Fig. 5.4, we observe that parametric model is essential to study the behavior of the COVID model. The dynamical behavior of our proposed COVID-19 model with respect time for different values of p is represented by Fig. 5.5 and Fig. 5.6, respectively. 0 50 100 150 200 250 300 350 400 450 500 0 1 2 3 4 5 6 7 x 10 7 Time A c ti v e i n fe c te d ( I S ) p=0.0 p=0.2 p=0.4 p=0.6 p=0.8 p=1.0 p=0.5 Fig. 5.5. Dynamical behavior with time for different values of p for the active infected population From Fig. 5.5, we observe that the pick of height occurs when p = 1 and which is nearly equal to 6× 107. It is also noticed that when p is gradually decreases form 1 then the number of height active individual is also decreases. For p = 0 the infective curve attains its pick level Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 83 0 100 200 300 400 500 2 4 6 8 10 12 x 10 8 Time S u s c e p ti b le ( S ) 0 100 200 300 400 500 0 1 2 3 4 5 x 10 7 Time E x p o s e d ( E ) 0 100 200 300 400 500 0 2 4 6 x 10 7 Time A s y m p to m a ti c ( I A ) 0 100 200 300 400 500 0 2 4 6 8 10 x 10 8 Time R e c o v e re d ( R ) p=0.0 p=0.2 p=0.4 p=0.6 p=0.8 p=1.0 p=0.5 i ii iviii Fig. 5.6. Dynamical behavior with time for exposed (E), susceptible (S), asymptomatic infected (IA), recovered (R) population curve for different values of p at nearly equal to 3× 107. Therefore, for our COVID-19 parametric model the total active number of infected individuals lies in the interval [3× 107, 6× 107]. Again, after reaching the height pick of the infection the active infective curve is gradually decreases with time for any values of p. Moreover, the curve of population E and IA shows the similar behavior as infective (IS) curve for different values of p with respect to time t. On the contrary, Susceptible (S) population curve shows an opposite behavior for different values of p. In this figure, we observe that as p decreases from 1 the susceptible curve gradually increases. Finally, we observed that as p increases from 0 the recovered population curve (R) is also increases, i.e., the number of recovered population gradually increased. Now, we solve the imprecise optimality system numerically and the obtained results are plotted graphically. The proposed optimality system is a two point boundary value system with divided boundary conditions at times t = 0 and t = tf . Here, we have solved this two- point boundary value optimality problem for tf = 25 in the time unit day at which the treatment is stopped. As in the objective functional (4.10), the different populations and control functions are different scales, we are balancing the objective functional (4.10) by opting the weight constants G1, G2 , G3, K1 and K2 given in Table 5.4. Table 5.4. Values of the weight constants weight constant Value G1 0.1 G2 0.2 G3 0.2 K1 300000 K2 300000 The numerical simulation of the control problem is performed using the parameter values and weight constants specified in Table 5.1 and Table 5.4, respectively. Here, we search for two optimal control functions u1 and u2, minimize the objective functional (J) provided by (4.10). The graph of optimal control for u1 and u2 is presented in Fig. 7 for different p. 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time u 1 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time u 2 p=0.0 p=1.0 p=0.5 p=0.0 p=1.0 p=0.5 iii Fig. 5.7. Optimal control of two controls, i.e., treatment control (u 1 and u 2 ) using the parameter values given in Table 3 with G1 = 0.1, G2 = 0.2, G3 = 0.2, K1 = 300000, K2 = 300000 for different values of p. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 84 P.K. SANTRA ET AL. From the optimal control diagrams (Fig. 5.7), one could conclude that the total effort should be given in both controls at the start of the disease. The system’s dynamical behavior under control and without control is presented in Fig. 5.8 for different p. 0 5 10 15 20 25 1.1985 1.199 1.1995 1.2 1.2005 1.201 x 10 9 Time P o p u la ti o n S 0 5 10 15 20 25 0 2 4 6 x 10 5 Time P o p u la ti o n E 0 5 10 15 20 25 0 2 4 6 x 10 5 Time P o p u la ti o n I A 0 5 10 15 20 25 0 2 4 6 x 10 5 Time P o p u la ti o n I S p=0.5 p=0.5 p=0.0 p=0.0 p=1.0 p=1.0 i iii iv ii Fig. 5.8. Time series plots of (i) susceptible population, (ii) exposed population, (iii) infective population in asymptomatic phase, (iv) infective population in symptomatic phase with two controls and without control using the parameter values given in Table 5.1 and the values of the weight constants given in Table 5.4. The solid lines of Fig. 5.8 are presented as the population with control, and the dotted lines are presented as the population without control. Fig. 5.8 shows that using control, exposed class (E), asymptomatic class (IA) and symptomatic class (IS) are decreasing. This study indicates that the two controls can control the disease successfully. In Fig. 5.9 and Fig. 5.10, we have presented the comparative study without control (u1 = u2 = 0) through the time series diagrams of every population class, also for only first treatment control (u1 ̸= 0, u2 = 0) and only second treatment control (u1 = 0, u2 ̸= 0) for different values of p. 0 5 10 15 20 25 1.1985 1.199 1.1995 1.2 1.2005 x 10 9 Time P o p u la ti o n S 0 5 10 15 20 25 1 2 3 4 5 6 x 10 5 Time P o p u la ti o n E 0 5 10 15 20 25 0 1 2 3 4 5 6 x 10 5 Time P o p u la ti o n I A 0 5 10 15 20 25 1 2 3 4 5 6 x 10 5 Time P o p u la ti o n I S p=0.5 p=0.5 p=0.0 p=0.0 p=1.0 p=1.0 iii i ii iv Fig. 5.9. Time series plots of (i) susceptible, (ii) exposed, (iii) infected in asymptomatic phase, (iv) infected in symptomatic phase populations with u1 control and without control using input values given in Table 5.1 and Table 5.4. 0 5 10 15 20 25 1.1985 1.199 1.1995 1.2 1.2005 x 10 9 Time P o p u la ti o n S 0 5 10 15 20 25 1 2 3 4 5 6 x 10 5 Time P o p u la ti o n E 0 5 10 15 20 25 1 2 3 4 5 6 x 10 5 Time P o p u la ti o n I A 0 5 10 15 20 25 0 2 4 6 x 10 5 Time P o p u la ti o n I S p=0.5 p=0.5 p=0.0 p=0.0 p=1.0 p=1.0 i ii iviii Fig. 5.10. Time series plots of (i) susceptible, (ii) exposed, (iii) asymptomatic phase of infected, (iv) symptomatic phase of infected populations with u2 control and without control using the parameter values given in Table 5.1 and the values of the weight constants given in Table 5.4. The solid lines of Fig. 5.9 are presented as the population with u1 control and the dotted lines are presented as the population without control. In Fig. 5.10, the solid lines are presented as the population with u2 control and dotted lines are presented as the population without Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 85 control. These figures clearly show that using both controls together is more effective to control the disease compared to without any control and single controls. From the above observations, it can conclude that both the treatment controls together yield a relatively noticed result on the controlling of the disease. 6. CONCLUSIONS This paper considered a pandemic disease in recent years, the so-called COVID-19 infection. This study developed a five-compartmental epidemic model and investigated the dynamical behavior of this model. Most of the COVID-19 infection models are generally based on the assumption that the model’s parameters are precisely known. However, the scenario is different since it is impossible to know all the parameter values specifically. In this paper, we developed a method to discuss the dynamical behavior of COVID-19 epidemic model with imprecise parameters by considering that the model’s coefficients are ambiguous for the lack of precise numerical information. By using the concept of the next-generation matrix method, we have found R0 = Λ1−p l Λp u µ1−p l µp u [α1−p l αp uσρ 1−p l ρpuC + β1−p l βp uρ 1−p l ρpu{σ(1− k)r1−p l rpu +B(1− σ)}] ABC as BRN of the system, which helps us to determine the dynamical behavior of the system for all p ∈ [0, 1]. The system (2.5) is locally asymptotically stable in the DFE E0 when R0 < 1 for all p ∈ [0, 1]. The endemic equilibrium E1 exists when R0 > 1 and p ∈ [0, 1] and the system becomes unstable at E0 but under some conditions locally asymptotically stable at E1. The most important part of this study is to establish an optimal control for the pandemic model to minimize the asymptomatic and symptomatic populations and also to minimize the cost of treatment under impreciseness. We have considered two types of treatment as two controls to reduce the spread of the disease. The first treatment control u1, and the second treatment control u2 are designed in such a way that they minimize the objective functional as given in equation (4.10). The main mathematical findings for the dynamical behavior of the COVID-19 pandemic are also numerically verified using MATLAB. A comparison between the dynamical behavior of every population with three different control options is presented graphically for different values of p in Fig. 5.7 to Fig. 5.10, which clearly shows the effectiveness of the treatment control. In addition, it is precisely shown that optimal control is much more effective in reducing the number of asymptomatic and symptomatic individuals, implying that treatment can successfully control the spread of the disease. Furthermore, it was noticed that both treatments are essential at the beginning of the outbreak to prevent the spread of the disease. Thus, our analytical and numerical studies show that the optimal control policy is a beneficial technique to reduce the spread of COVID-19 infection. The proposed pandemic model is based on the effects of COVID-19 infection in a population. In most of the COVID-19 mathematical models, the model parameters are taken as constants. However, in reality, it is not true, as the parameters depend on the environmental conditions; therefore, they are not constant. To avoid these difficulties, in our present model we consider that the parameters of the models are imprecise and present them by interval number. Then, using the parametric functional form of the interval number, we study the dynamical behavior of the imprecise model. It is also true that a time lag between susceptible individuals may be infected, hence, the present model can be extended by incorporating the time lag into the system, which is kept for future work consideration. Controlling COVID-19 spread is now a challenging and significant issue. Therefore, predicting and identifying appropriate strategies and minimum cost prevention programs to stop the spread of the virus are the government’s primary goals. This study is a step toward identifying the parameters of interest for future studies. However, a tremendous amount of Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 86 P.K. SANTRA ET AL. contribution is needed to inform and assist researchers and policy makers in focusing on prevention and available treatment resources for maximum effectiveness. REFERENCES 1. Kumar, R.P., Santra, P.K. & Mahapatra, G.S. (2023). Global stability and analysing the sensitivity of parameters of a multiple-susceptible population model of SARS-CoV-2 emphasising vaccination drive, Math. Comput. Simul., 203, 741–766. 2. Santra, P.K., Ghosh, D., Mahapatra, G.S. & Bonyah, E. (2022). Mathematical Analysis of Two Waves of COVID-19 Disease with Impact of Vaccination as Optimal Control, Comput. Math. Methods. Med., 2022, 2684055. 3. Kim, S., Byun, J.H. & Jung, H. (2019). Global stability of an SEIR epidemic model where empirical distribution of incubation period is approximated by Coxian distribution, Adv. Differ. Equ., 469. 4. Kumar, R.P., Basu, S., Santra, P.K., Ghosh, D. & Mahapatra, G.S. (2022). Optimal control design incorporating vaccination and treatment on six compartment pandemic dynamical system, Results in Control and Optimization, 7, 100115. 5. Basu, S., Kumar, R.P., Santra, P.K., Mahapatra, G.S. & Elsadany, A.A. (2022). Preventive control strategy on second wave of COVID-19 pandemic model incorporating lock-down effect, Alex. Eng. J., 61(9), 7265–7276. 6. Kumar, R.P., Basu, S., Ghosh, D., Santra, P.K., & Mahapatra, G.S. (2022). Dynamical analysis of novel COVID-19 epidemic model with non-monotonic incidence function, J. Public Aff., 22, e2754. 7. Cohen, J., & Normile, D. (2020). New SARS-like virus in China triggers alarm, Science, 367(6475), 234–235. 8. World Health Organization. (2020, May 8). Naming the coronavirus disease (COVID-19) and the virus that causes it. [Online]. Available: https://www.who.int/emergencies/diseases/novel-coronavirus-2019/technical- guidance/naming-the-coronavirus-disease-(COVID-2019)-and-the-virus-that-causes-it. 9. World Health Organization. (2020, May 8). Coronavirus dis- ease 2019 (COVID-19), situation report 109. [Online]. Available: https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports. 10. Pai, C., Bhaskar, A. & Rawoot, V. (2020). Investigating the dynamics of COVID-19 pandemic in India under lockdown, Chaos Soliton Fract., 138, 109988. 11. Djilali, S. & Ghanbari, B. (2020). Coronavirus pandemic: A predictive analysis of the peak outbreak epidemic in South Africa, Turkey, and Brazi, Chaos Soliton Fract., 138, 109971. 12. Ayinde, K., Lukman, A.F., Rauf, R.I., Alabi, O.O., Okon, C.E. & Ayinde, O.E. (2020). Modeling Nigerian COVID-19 cases: A comparative analysis of models and estimators, Chaos Soliton Fract., 138, 109911. 13. Kırbaş, I., Sözen, A., Tuncer, A.D. & Ş.Kazancıoğlu, F. (2020). Comparative analysis and forecasting of COVID-19 cases in various European countries with ARIMA, NARNN and LSTM approaches, Chaos Soliton Fract., 138, 110015. 14. Păcurar, C.M. & Necula, B.R. (2020). An analysis of COVID-19 spread based on fractal interpolation and fractal dimension, Chaos Soliton Fract., 139, 110073. 15. COVID-10India, (2020, May 8). India COVID-19 tracker. [Online]. Available: https://www.COVID19india.org/. 16. Pan, Q., Gao, T. & He, M. (2020). Influence of isolation measures for patients with mild symptoms on the spread of COVID-19, Chaos Soliton Fract., 139, 110022. 17. Panwar, H., Gupta, P.K., Siddiqui, M.K., Morales-Menendez, R. & Singh, V. (2020). Application of deep learning for fast detection of COVID-19 in X-Rays using nCOVnet, Chaos Soliton Fract., 138, 109944. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 87 18. Huang, Y., Wu, Y. & Zhang, W. (2020). Comprehensive identification and isolation policies have effectively suppressed the spread of COVID-19, Chaos Soliton Fract., 139, 110041. 19. Silva, P.C.L., Batista, P.V.C., Lima, H.S., Alves, M.A., Guimarães, F.G. & Silva, R.C.P. (2020). An agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions, Chaos Soliton Fract., 139, 110088. 20. Ferguson, N., Laydon, D., Gilani, G.N., Imai, N., Ainslie, K., et al. (2020, March 16). Report 9: Impact of non-pharmaceutical interventions (npis) to reduce COVID19 mortality and healthcare demand. [Online]. Available: https://www.imperial.ac.uk/media/imperial-college/medicine/mrc-gida/2020-03- 16-COVID19-Report-9.pdf. 21. Marquioni, V.M. & Aguiar, M.A.M.D., (2020). Quantifying the effects of quarantine using an IBM SEIR model on scalefree networks, Chaos Soliton Fract., 138, 109999. 22. Fredj, H.B. & Chérif, F. (2020). Novel Corona virus disease infection in Tunisia: Mathematical model and the impact of the quarantine strategy, Chaos Soliton Fract., 138, 109969. 23. Tang, B., Bragazzi, N.L., Li, Q., Tang, S., Xiao, Y. & Wu, J. (2020). An updated estimation of the risk of transmission of the novel coronavirus (2019-ncov), Infectious Disease Modelling, 5, 248–255. 24. Anguelov, R., Banasiak, J., Bright, C., Lubuma, J. & Ouifki, R. (2020). The big unknown: The asymptomatic spread of COVID-19, Biomath., 9, 2005103. 25. Okuonghae, D. & Omame, A. (2020). Analysis of a mathematical model for COVID-19 population dynamics in Lagos, Nigeria, Chaos Soliton Fract., 139, 110032. 26. Tang, B., Wang, X., Li, Q., Bragazzi, N.L., Tang, S., Xiao, Y. & Wu, J. (2020). Estimation of the transmission risk of the 2019-ncov and its implication for public health interventions, Journal of Clinical Medicine, 9(2), 462. 27. Askar, S.S., Ghosh, D., Santra, P.K., Elsadany, A.A. & Mahapatra, G.S. (2021). A fractional order SITR mathematical model for forecasting of transmission of COVID-19 of India with lockdown effect, Results in Physics, 24, 104067. 28. Cadoni, M. & Gaeta, G. (2020). Size and timescale of epidemics in the SIR framework, Physica D, 411, 132626. 29. Zhang, Y., Yu, X., Sun, H., Tick, G.R., Wei, W., & Jin, B. (2020). Applicability of time fractional derivative models for simulating the dynamics and mitigation scenarios of COVID-19, Chaos Soliton Fract., 138, 109959. 30. Higazy, M. (2020). Novel fractional order SIDARTHE mathematical model of COVID- 19 pandemic, Chaos Soliton Fract., 138, 110007. 31. Wu, J.T., Leung, K. & Leung, G.M. (2020). Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in Wuhan, China: a modelling study, The Lancet, 395(10225), 689–697. 32. Read, J.M., Bridgen, J.R., Cummings, D.A., A. Ho, A. & Jewell, C.P. (2020). Novel coronavirus 2019-ncov: early estimation of epidemiological parameters and epidemic predictions, Philosophical Transactions of the Royal Society B: Biological Sciences, 376(1829), 1–9. 33. Pal, D., Ghosh, D., Santra, P.K. & G.S. Mahapatra, (2022). Mathematical analysis of a COVID-19 epidemic model by using data driven epidemiological parameters of diseases spread in India, Biophysics, 67, 231–244. 34. Nandi, S.K., Jana, S., Manadal, M. & Kar, T.K. (2018). Analysis of a fuzzy epidemic model with saturated treatment and disease transmission, Int J Biomath, 11(1), 1850003- 1–1850003-20. 35. Cai, Y., Kang, Y. & Wang, W. (2017). A stochastic SIRS epidemic model with nonlinear incidence rate, Appl. Math. Comput., 305, 221–40. 36. Lahrouz, A., Omari, L. & Kiouach, D. (2011). Global analysis of a deterministic and stochastic nonlinear SIRS epidemic model, Nonlinear Anal Model Control, 16, 59–76. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 88 P.K. SANTRA ET AL. 37. Driessche, P.V.D. & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180, 29–48. 38. Çakan, S. (2020). Dynamic analysis of a mathematical model with health care capacity for COVID-19 pandemic, Chaos Soliton Fract., 139, 110033. 39. Ng, K.Y. & Gui, M.M. (2020). COVID-19: Development of a robust mathematical model and simulation package with consideration for ageing population and time delay for control action and resusceptibility, Physica D, 411, 132599. 40. Al-Asuoad, N., Alaswad, S., Rong, L. & Shillor, M. (2016). Mathematical model and simulations of MERS outbreak: Predictions and implications for control measures, Biomath, 5, 1612141. 41. Zhao, H. & Feng, Z. (2020). Staggered release policies for COVID-19 control: Costs and benefits of relaxing restrictions by age and risk, Chaos Soliton Fract., 138, 108405. 42. Abbasi, Z., Zamani, I., Mehra, A.H.A, Shafieirad, M. & Ibeas, A. (2020). Optimal Control Design of Impulsive SQEIAR Epidemic Models with Application to COVID- 19, Chaos Soliton Fract., 139, 110054. 43. Mbazumutima, V., Thron, C. & Todjihounde, L. (2019). Enumerical solution for optimal control using treatment and vaccination for an SIS epidemic model, Biomath, 8, 1912137. 44. Mahapatra, G.S. & Mandal, T.K. (2012). Posynomial Parametric Geometric Programming with Interval Valued Coefficient, Journal of Optimization Theory and Applications, 154, 120–132. 45. Santra, P.K. & Mahapatra, G.S. (2020). Dynamical study of of Discrete-time Prey- predator Model with Constant Prey Refuge under Imprecise Biological Parameters, Journal of Biological Systems, 28(3), 681–699. 46. Brauer, F. & Castillo-Chavez, C. (2001). Mathematical Models in Population Biology and Epidemiology. Berlin, Germany: Springer. 47. Driessche, P.V.D. & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180, 29–48. 48. Castillo-Chavez, C., Feng, Z. & Huang, W. (2002). On the computation of R0 and its role on Global Stability. In: C. Castillo-Chavez, S. Blower, P. Driessche, D. Kirschner, & A.-A. Yakubu (Eds.), Mathematical approaches for emerging and reemerging infectious diseases: an introduction. 49. Kot, M. (2001). Elements of Mathematical Ecology. Cambridge, UK: Cambridge University Press. 50. Castillo-Chevez, C., & Feng, Z. (1998). Global stability of an agestructure model for TB and its applications to optimal vaccination strategies, Math Biosci, 151, 135–154. 51. Fleming, D.T. & Wasserheit, J.N. (1999). From epidemiological synergy to public health policy and practice: the contribution of other sexually transmitted diseases to sexual transmission of HIV infection, Sex Transm. Infect., 75(4), 3–17. 52. Lukes, D.L. (1982). Differential Equations: Classical to Controlled, Mathematics in Science and Engineering. New York, NY: Academic Press. 53. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V. & Mishchenko, E.F. (1986). The mathematical theory of optimal processes. London, UK: Gordon and Breach Science. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) COVID-19 SPREAD MODELING UNDER UNCERTAINTY 89 Appendix A: The proof of theorem-6 The variational matrix of system (2.5) at DFE E0 is given by ME0 =  −µ1−p l µp u 0 −α1−p l αp uΛ 1−p l Λp u µ1−p l µu)p −β1−p l βp uΛ 1−p l Λp u µ1−p l µp u 0 0 −A α1−p l αp uΛ 1−p l Λp u µ1−p l µp u β1−p l βp uΛ 1−p l Λp u µ1−p l µp u 0 0 σρ1−p l ρpu −B 0 0 0 (1− σ)ρ1−p l ρpu (1− k)r1−p l rpu −C 0 0 0 kr1−p l rpu m1−p l mp u −µ1−p l µp u  Therefore, eigenvalues of the characteristic equation of ME0 are −µ1−p l µp u, −µ1−p l µp u and the solution of the cubic equation, P (λ) ≡ λ3 + A1λ 2 + A2λ+ A3 = 0 (A.1) where A1 = A+B + C, A2 = AB(1−R0) + AC(1−R0) +BC + Λ1−p l Λp u µ1−p l µp u[ ρ1−p l ρpuβ 1−p l βp u (1− σ)B + (1− k)σr1−p l rpu C + σρ1−p l ρpu α1−p l αp uC + (1− k)β1−p l βp ur 1−p l rpu B ] and A3 = ABC(1−R0) Now, it is easily noted that A1 > 0, A2 > 0, A3 > 0 if R0 < 1. After some simplifications, we get A1A2 −A3 = (A2B +AB2 +A2C +AC2 +ABC)(1−R0) + (A+B + C)( BC + Λ1−p l Λp uρ 1−p l ρpu µ1−p l µp u [ β1−p l βp u (1− σ)B + σ(1− k)r1−p l rpu C + σ α1−p l αp uC + (1− k)β1−p l βp ur 1−p l rpu B ]) Here, we can notice that, if R0 < 1 then A1A2 − A3 > 0. Therefore, by the Routh–Hurwitz Routh–Hurwitz criterion [49] it follows that P (λ) = 0 has negative real roots if R0 < 1.This completes the proof. Appendix B: The proof of theorem-8 The variational matrix of system (2.5) at E1(S ∗, E∗, I∗A, I ∗ S, R ∗) is given by, ME1 =  b11 0 b13 b14 0 b21 b22 b23 b24 0 0 b32 b33 0 0 0 b42 b43 b44 0 0 0 b53 b54 b55  where, b11 = − ( α1−p l αp uI ∗ A + β1−p l βp uI ∗ S ) − µ1−p l µp u, b13 = −α1−p l αp uS ∗, b14 = −β1−p l βp uS ∗, b21 = ( α1−p l αp uI ∗ A + β1−p l βp uI ∗ S ) , b22 = −A, b23 = α1−p l αp uS ∗, b24 = β1−p l βp uS ∗, b32 = Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 90 P.K. SANTRA ET AL. σρ1−p l ρpu, b33 = −B, b42 = (1− σ)ρ1−p l ρpu, b43 = (1− k)r1−p l rpu, b44 = −C, b53 = kr1−p l rpu, b54 = m1−p l mp u, b55 = −µ1−p l µp u. Therefore, the eigenvalues of the characteristic equation of ME1 are −µ1−p l µp u and the solution of the equation, Q (λ) ≡ λ4 +B1λ 3 +B2λ 2 +B3λ+B4 = 0 (B.1) where B1 = −(b11 + b22 + b33 + b44), B2 = b11b22 + b11b33 + b11b44 + b22b44 + b33b44 + b22b33 − b23b32 − b24b42, B3 = −b11b22b44 − b11b33b44 − b11b22b33 − b22b33b44 − b24b32b43 − b21b13b32 − b21b14b42 + b11b23b32 + b11b24b42 + b23b32b44 + b24b42b33, and B4 = b11b22b33b44 + b11b24b32b43 + b21b13b32b44 + b21b14b33b42 − b11b23b32b44 − b11b24b33b42 − b14b21b32b43. By the Routh–Hurwitz criterion [49] it follows that Q (λ) = 0 has negative real roots if Bi > 0 (i = 1, 2, 3, 4), D1 = B1 > 0, D2 = ∣∣∣∣B1 B3 1 B2 ∣∣∣∣ = B1B2 −B3 > 0, D3 =∣∣∣∣∣B1 B3 0 1 B2 B4 0 B1 B3 ∣∣∣∣∣ = B1B2B3 −B2 1B4 −B2 3 > 0. Therefore, the system (2.5) shows local asymptotic stability at E1 when R0 > 1, Bi > 0(i = 1, 2, 3, 4), B1B2 −B3 > 0 and B1B2B3 −B2 1B4 −B2 3 > 0. This completes the proof. Appendix C: Arithmetic operations on interval numbers using the concept of interval-valued functions [45] between A = [a, a] and B = [ b, b ] for p ∈ [0, 1] are as follows: Addition: A+B = [a, a] + [ b, b ] = [ a+ b, a+ b ] , then the interval-valued function for A+B is h (p) = c (1−p) L cpU where cL = a+ b > 0 and cU = a+ b. Subtraction: A−B = [a, a]− [ b, b ] = [ a− b, a− b ] , then the interval-valued function for A−B is h (p) = c (1−p) L cpU where cL = a− b > 0 and cU = a− b. Scalar Multiplication: εA = ε [a, a] = { [εa, εa] if ε ≥ 0 [αa, εa] if ε < 0 , provided a > 0 , then the interval-valued function of εA is h (p) = c (1−p) L cpU if ε ≥ 0 and h (p) = −d (1−p) L dpU if ε < 0, where c L = εa, c U = εa, dL = |ε|a and dU = |ε|a. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) Introduction Materials & Methods Variables and parameters involved in the Model Some Useful Definitions COVID-19 pandemic model formulation Crisp model Imprecise model with interval coefficient Theoretical study of the model Positivity of the solutions Invariant region Existence of Equilibria and Basic reproduction number Stability Behavior Optimal control for COVID-19 model Existence of an optimal control Characterization of the optimal control Numerical Simulations Conclusions