1 Introduction

First identified in Wuhan City, China, the COVID-19 disease (caused by the SARS-COV-2 Virus) is currently one of the worst pandemics in human history. The rapid spread of the virus throughout the world has caused severe economic, demographic and societal problems. During the first 12 months of its appearance, and given the lack of information on the virus, health authorities have opted for proactive measures to reduce the extent of the damage caused by this pandemic. Non-pharmaceutical measures, such as lockdowns, social distancing, quarantine and the wearing of masks (Perra 2021), were the main weapons available against the disease. Additionally, some treatments combining several drugs have been proposed to treat confirmed patients [(see Rismanbaf (2020)].

Despite these control measures, the reported infection and death cases have been stepped up across several countries. According to the World Health Organization (WHO), there have been 158,651,638 confirmed cases of COVID-19, including 3,299,764 deaths around the world, as of May 11, 2021 (Organization et al. 2020). In this context, the use of vaccines arises as an important alternative that can protect against COVID-19 and save the lives of millions. In early 2021, several companies began distributing the first COVID-19 vaccines. However, due to high demand and limited supply of vaccines, many countries are finding it very difficult to procure the necessary quantities and carry out their vaccination programmes.

It is well known that mathematical modelling can predict the evolution of a pandemic and help to put in place control strategies that take into consideration certain constraints related to the nature of the disease and the resources available to public authorities [see, e.g. Kumar and Srivastava (2017); Elhia et al. (2021); Zamir et al. (2021) and references therein]. To our knowledge, there is a scare on modelling studies taking into consideration the limitation on vaccine resources. The most popular work is done by Neilan and Lenhart (2010), which control the propagation of the disease by a limited number of vaccines over a period of time T while minimizing the number of infected people and a cost associated with vaccination. Based on this approach, Biswas et al. (2014) concentrated on the case where the supply of vaccines is limited at each unit of time, which can be modelled as a mixed state-control constraint. Recently, such approaches are also applied to diseases like Ebola and Influenza, see, e.g. Area et al. (2018), Wang et al. (2019), Zhang et al. (2020).

In the case of COVID-19, several studies have been proposed to understand the dynamics of this infectious disease and investigate the impact of some control measures. For example, the work done by Libotte et al. (2020) adopted the SIR model to describe the dynamic behaviour of the COVID-19 epidemic in China and proposed two optimal control problems (mono- and multi-objective) to determine a strategy for vaccine administration in COVID-19 pandemic treatment. Moore et al. (2021) used a mathematical model structured by age and UK region to investigate several vaccination scenarios and predict the possible long-term dynamics of SARS-CoV-2 when the UK vaccination programme is combined with multiple potential future relaxations of non-pharmaceutical interventions. Elhia et al. (2021) formulated a free terminal time optimal control problem that aims to find the optimal duration of a vaccination campaign using data from Morocco.

It must be emphasised that the models that existed in the literature differ on the choice of the compartmental model, on constraints enforced or on the functional cost. These models are generally governed by ordinary (Emary et al. 2021), partial (Miyaoka et al. 2019) or stochastic differential equations (Britton et al. 2019), and fractional-order derivatives [see, e.g. Sweilam et al. (2019), Boujallal (2021)]. In addition, the control problems formulated for COVID-19 have been solved using direct (Kohler et al. 2020) and indirect approaches (Grigorieva et al. 2021; Khajji et al. 2021) from control theory. For a detailed literature review of the different types of models proposed for COVID-19, the interested readers can consult (Perra 2021).

In this work, we focus on the control of the COVID-19 pandemic by setting and resolving a more realistic control problem that takes into account the limited vaccines supply and treatment. The novelty of our work consists of

  • Considering a mixed constraints control problem by adding the limited vaccines supply as constraint. To the best of our knowledge, our work is the first that takes the limitation of the COVID-19 vaccines resources as a constraint when formulating the control problem.

  • Proposing an extension of our previous work (Boujallal et al. 2021b) to the case with mixed constrains.

  • Using the set-valued approach which is predicated on viability theory (Aubin 2009) and set-valued analysis (Aubin and Frankowska 2009), we provide feedback controls via selections of a suitable constructed multifunction, which take into consideration constraints on the control variables and the vaccine limitation resources, and can reduce to zero asymptotically the number of infected people.

We recall that viability theory offers a solid framework of dynamical systems that obey constraints in a set-valued form. The main concern of viability theory is the computation of the set of initial states (known as a viability kernel) from which the evolution of a given system remains within a viability constraint set. Motivated by the difficulty of kernel computation, we use a set-valued approach that allows the characterization of the viability set in terms of a contingent cone (see definition 14). This approach is used for systems in finite and infinite dimensions [see Boujallal and Kassara (2015), Boujallal and Kassara (2020)]. More recently, this approach is employed in (Boujallal et al. 2021a, b) to deal with a control problem of a general class of epidemiological models where the objective is to eradicate the infectious from a population, while Elhia et al. (2020) treated this problem for the SEIR model of COVID-19 with saturated treatment function and screening effort.

The remainder of the paper is organised as follows: In Sect. 2, we present the SEIR COVID-19 model and set the control problem we will investigate. Section 3 will be concerned with presenting the set-valued approach, and Sect. 4 to characterise a continuous selection that provides us with control strategies to deal with the COVID-19 pandemic. To show the efficiency of our approach, we propose in Sect. 5 numerical simulations and conclude this paper in Sect. 6.

2 Mathematical Model and Control Problem

2.1 Mathematical Model

The COVID-19 model under consideration in the present paper describes the dynamics of a population divided into four different compartments

  • Susceptible individuals (S) are those vulnerable to become infected;

  • Exposed individuals (E), who have been infected but are not infectious and do not show symptoms of the disease;

  • Infectious individuals (I), who may transmit the infection;

  • Recovered or Immune individuals (R).

Thus, the total population at time t is presented by

$$\begin{aligned} N(t)=S(t)+E(t)+I(t)+R(t). \end{aligned}$$

This model is formulated based on the following parameters

  • All recruitments are into the susceptible compartment and occur at a constant rate \(\Lambda\);

  • d denotes the natural death rate;

  • The transmission of the disease occurs following an adequate contact between a susceptible and infectious. Due to the nonlinear contact dynamics in the population, we use the incidence function \(c\,\dfrac{SI}{N}\) to indicate successful transmission of the disease, where c denotes the effective contact rate with infectious individuals in compartment I;

  • The exposed individuals become infectious at a rate e;

  • g is the rate at which infectious individuals recover;

  • a denotes the death rate due to the disease.

Here, all the parameters are positive constants.

The dynamics of the model is governed by the following nonlinear system

$$\begin{aligned} \left\{ \begin{array}{lll} \dot{S}{}\,=\,{}\Lambda -c\,\dfrac{SI}{N}-\mathrm{d}S, \\ \dot{E}{}\,=\,{} c\dfrac{SI}{N}-(e+d)E,\\ \dot{I}{}\,=\,{} eE-(g+a+d)I,\\ \dot{R}{}\,=\,{}gI-dR, \end{array} \right. \end{aligned}$$
(1)

with initial conditions

$$\begin{aligned} \begin{array}{ll} S(0)=S_0\ge 0,~ E(0)=E_0\ge 0,\\ I(0)=I_0\ge 0,~ R(0)=R_0\ge 0. \end{array} \end{aligned}$$
(2)

2.2 Control Problem Under Mixed Constraints

In this section, we introduce two control strategies to eradicate the COVID-19 infectious from a specified population, namely vaccination and treatment strategies. Motivated by realistic condition, the limitation of vaccine resources at each time instant is also considered in this paper. Thus, this problem is solved in the framework of set-valued approach as stated in Sect. 3, and a modified dynamics of the SEIR model (1) with vaccination and treatment is proposed as follows:

$$\begin{aligned} \left\{ \begin{array}{lll} \dot{S}{}\,=\,{}\Lambda -c\,\dfrac{SI}{N}-\mathrm{d}S-u_1S, \\ \dot{E}{}\,=\,{} c\dfrac{SI}{N}-(e+d)E,\\ \dot{I}\,=\,{} eE-(g+a+d)I-u_2I,\\ \dot{R}{}\,=\,{}gI-dR+u_1S+u_2I, \end{array} \right. \end{aligned}$$
(3)

subject to control constraints

$$\begin{aligned} 0\le u_i \le u_i^{\max } \le 1~ for~ i=1,2, \end{aligned}$$
(4)

and mixed state-control constraints

$$\begin{aligned} u_1S\le V_0, \end{aligned}$$
(5)

with \(V_0\) is an upper bound taking positive values.

The control in system (3) is represented by \(u=(u_1,u_2)^\prime\), where \(u_1\) stands for the fraction of vaccinated individuals and \(u_2\) represents the fraction of individuals that is put under treatment. The control u must range in the subset,

$$\begin{aligned} \mathcal {U}\doteq \prod _{i=1}^{2}[0, u_i^\mathrm{max}],~~\hbox {where}~u_i^\mathrm{max}\le 1. \end{aligned}$$

Therefore, the control problem we have to deal with is stated as follows:

For each \((S_0, E_0, I_0, R_0)\in \Gamma \subset {\mathbb {R}}^4_+\), find a control \(\bar{u}\) such that:

$$\begin{aligned} \bar{u}: [0,\infty )\rightarrow \mathcal {U}, \end{aligned}$$
(6a)
$$\begin{aligned} \bar{u}_1(t)\bar{S}(t)\le V_0,~~\hbox {for all}~t\ge 0, \end{aligned}$$
(6b)

and,

$$\begin{aligned} \displaystyle {\lim _{t\rightarrow \infty }\bar{I}(t)=0}, \end{aligned}$$
(6c)

where \((\bar{S},\bar{E},\bar{I},\bar{R})\) denotes a solution of system (3) for the control \(\bar{u}\).

3 The Set-Valued Approach

Here, we are going to provide a brief outline of the set-valued control approach mentioned above, see Boujallal (2022) for more mathematical details. It involves the following class of nonlinear augmented control system,

$$\begin{aligned}&\displaystyle {\dot{x}=f(x,y,u),}~~(x,y,u)\in \mathcal {K}, \end{aligned}$$
(7a)
$$\begin{aligned}&\displaystyle {\dot{y}=h(x,y,u),} \end{aligned}$$
(7b)
$$\begin{aligned}&\displaystyle {\dot{u}=v-\delta u,} \end{aligned}$$
(7c)

where f and h represent sufficiently smooth functions from \({\mathbb {R}}^n\times {\mathbb {R}}^m\times {\mathbb {R}}^p\) to \({\mathbb {R}}^n,\) and \({\mathbb {R}}^m,\) respectively, for integers n, m and p. The state (x, y) evolves in the subset \(S\doteq {\mathbb {R}}^n_+\times {\mathbb {R}}^m_+\), the control u takes values within constraint subset, \(\mathcal {U}\doteq \prod _{i=1}^{p}[0,u_i^\mathrm{max}]\), where \(u_i^\mathrm{max}\) are positive numbers, and \(\mathcal {K}\) stands for a mixed state-control constraint set. Here, v stands for the control and has values in closed subset

$$\begin{aligned} \displaystyle \mathcal {V}=\prod _{i=1}^{i=2}[0, \delta u_i^{\max }]. \end{aligned}$$
(8)

In the sequel, we need to assume the following assumptions:

Assumption 1

$$\begin{aligned} ~\begin{array}{|llll} x_i=0 \;\hbox {and}\; y\ge 0\implies f_i(x,y,u)\ge 0,\;\hbox {for all}\;\;u\in \mathcal {U}\\ y_j=0 \;\hbox {and}\; x\ge 0\implies h_j(x,y,u)\ge 0,\;\;\hbox {for all}\;\;u\in \mathcal {U}, \end{array} \end{aligned}$$
(9)

for  \(i=1,...,n\)  and  \(j=1,...,m.\)

Assumption 2

  There exist constants \(c_1, c_2>0\) such that, for all (x, y, u) in \(\mathcal {D}_\alpha\),

$$\begin{aligned} \Vert f(x,y,u)\Vert \le c_1(\Vert x\Vert +\Vert y\Vert +\Vert u\Vert +1), \end{aligned}$$

and

$$\begin{aligned} \Vert h(x,y,u)\Vert \le c_2(\Vert x\Vert +\Vert y\Vert +\Vert u\Vert +1). \end{aligned}$$

Remark 1

It is noteworthy to mention that, under Assumption 1, all trajectories of system (), issued from \((x_0,y_0)\in {\mathbb {R}}^n_+\times {\mathbb {R}}^m_+\), are positively invariant; i.e. \((\bar{x}(t),\bar{y}(t))\in {\mathbb {R}}^n_+\times {\mathbb {R}}^m_+\) for all \(t\ge 0\), while assumption 2 guarantees the existence of global solution, i.e. on \([0,\infty )\).

Remark 2

The solution \(\bar{u}\) of (7c) related to control \(v\in \mathcal {V}\) satisfies the inequality:

$$\begin{aligned} u_j^0\exp (-\alpha t)\le \bar{u}_j(t)\le u_j^{\max }-(u_j^{\max }-u_j^0)\exp (-\alpha t), \end{aligned}$$

for \(j=0,\dots ,p\). Then, for each \(u_0\in \mathcal {U}\),

$$\begin{aligned} \bar{u}(t)\in \mathcal {U}~ \hbox {for all} ~t\ge 0. \end{aligned}$$
(10)

Definition 1

We call a control strategy any function \(\bar{v}:[0,\infty )\rightarrow \mathcal {V}\) such that system () has a global solution \((\bar{x},\bar{y},\bar{u})\) with values in \(\mathcal {K}\cap {\mathbb {R}}^4_+\times \mathcal {U}\), and satisfying \(\displaystyle \lim _{t\rightarrow \infty }\bar{y}(t)=0\).

Definition 2

(Boujallal 2022) Let \(\Omega\) stand for a non empty subset of \({\mathbb {R}}^m\) and, \(\varphi : {{\mathbb {R}}^m}\times {{\mathbb{R}}^m}\rightarrow {\mathbb {R}},\) a \(\mathcal {C}_1\) real-valued function. It is called an \(\Omega\)-Lyapunov function if the following holds,

$$\begin{aligned} \begin{array}{llll} s : [0~\infty ) \rightarrow \Omega , \hbox { differentiable, and }\\ \varphi (s(t),\dot{s}(t))\le 0 \hbox { for all } t\ge 0, \end{array} \implies \lim _{t\rightarrow \infty } s(t)= 0. \end{aligned}$$
(11)

Inspired by Lyapunov functions given in Definition 2, we introduce the subset

$$\begin{aligned} \displaystyle \mathcal {D}_\varphi \doteq \left\{ (x,y,u)\in \mathcal {K} ~\mid ~ \varphi (y,h(x,y,u))\le 0 \right\} , \end{aligned}$$
(12)

and its related feedback map given by

$$\begin{aligned} \displaystyle \mathcal {G}_\varphi (x,y,u)~\doteq \left\{ v\in \mathcal {V}~\mid ~\right.&\nonumber \\ \displaystyle \left. (f(x,y,u), h(x,y,u), v-\delta u)\in T_{\mathcal {D}_\varphi }(x,y,u)\right\} , \end{aligned}$$
(13)

where

$$\begin{aligned} \displaystyle {T_{\mathcal {D}_\varphi }(z)=\{q~~ \mid \liminf _{\varepsilon \downarrow 0}\frac{d(z+\varepsilon q, {\mathcal {D}_\varphi })}{\varepsilon }=0\},} \end{aligned}$$
(14)

stands for the contingent cone to subset \(\mathcal {D}_\varphi\) at point \(z\in \mathcal {D}_\varphi\). Therefore, feedback control strategies can be given as continuous selections \(v=\zeta (x,y,u)\) of the feedback map \(\mathcal {G}_\varphi\) given by (13). Recall that \(\zeta\) is a selection of \(\mathcal {G}_\varphi\) whenever \(\zeta (x,y,u)\in \mathcal {G}_\varphi (x,y,u)\) for all \((x,y,u)\in \mathcal {D}_\varphi\).

According to Boujallal (2022), the feedback map given in (13) can be expressed for each \((x,y,u)\in \mathcal {D}_\varphi\), as follows:

$${\mathcal{G}}_{\varphi } (x,y,u) = ~\left| {\begin{array}{*{20}l} {{\mathcal{G}}(x,y,u)} & {{\text{if}}\,{\text{ }}\varphi (y,h(x,y,u)) < 0,} \\ {{\mathcal{C}}_{\varphi } (x,y,u)} & {{\text{if }}\,\varphi (y,h(x,y,u)) = 0,} \\ \end{array} } \right.$$
(15)

where

$$\begin{aligned} \mathcal {G}(x,y,u)\doteq \left\{ v\in \mathcal {V}~\mid ~ \right. \nonumber \\ \left. (f(x,y,u),h(x,y,u),v-\delta u)\in T_{\mathcal {K}}(x,y,u)\right\} , \end{aligned}$$
(16)

and

$$\begin{aligned} \mathcal {C}(x,y,u)\doteq \left\{ v\in \mathcal {G}(x,y,u)~\mid ~\right. \nonumber \\ \left. \ell _\varphi (x,y,u)+\left\langle m_\varphi (x,y,u),v\right\rangle \le 0 \right\} , \end{aligned}$$
(17)

for all \((x,y,u)\in \mathcal {K}.\) Here, the functions \(\ell _\varphi\) and \(m_\varphi\) are given by

$$\begin{aligned}& \ell _\varphi (x,y,u)\doteq \left\langle \nabla _z\varphi (y,h(x,y,u)),\nabla _xh(x,y,u)f\right\rangle\\ &\quad+\left\langle \nabla _z\varphi (y,h(x,y,u)),\nabla _yh(x,y,u)h\right\rangle \\& \quad-\delta \left\langle \nabla _z\varphi (y,h(x,y,u)),\nabla _uh(x,y,u) u\right\rangle \\ &\quad+\left\langle \nabla _y\varphi (y,h(x,y,u)),h(x,y,u)\right\rangle, \end{aligned}$$
(18)
$$\begin{aligned} m_\varphi (x,y,u)\doteq (\nabla _uh(x,y,u))^\prime \nabla _z\varphi (y,h(x,y,u)). \end{aligned}$$
(19)

Assumption 3

For all \((x,y,u)\in \mathcal {D}_\varphi\), there exists \(v \in \mathcal {G}(x,y,u)\) which satisfies the following statement

$$\begin{aligned} \ell _\varphi (x,y,u)+\left\langle v, m_\varphi (x,y,u)\right\rangle <0, \end{aligned}$$
(20)

where functions \(\ell _\varphi\) and \(m_\varphi\) are, respectively, given by Eqs. (18) and (19).

Remark 3

It should be also noted that assumption 3 plays a central rule in proving the existence of continuous selections of the feedback map \(\mathcal {G}_\varphi\).

Then, we are in a statement to prove the following main results.

Theorem 1

Assume that both assumptions 2 and 3 are satisfied. Then, any continuous selection of the feedback map \(\mathcal {G}_\varphi\) provides feedback control strategies for each \((x_0,y_0,u_0)\in \mathcal {D}_\varphi \cap {\mathbb {R}}^4_+\times \mathcal {U}\).

Proof

Let \(\zeta :\mathcal {D}_\varphi \rightarrow \mathcal {V}\) be a continuous selection of \(\mathcal {G}_\varphi\), then subset \(\mathcal {D}_\varphi\) is viable under system ().

This means that, for each \((x_0, y_0,u_0)\in \mathcal {D}_\varphi\), the system () has a solution \((\bar{x},\bar{y},\bar{u}):[0,\infty )\rightarrow \mathcal {D}_\varphi ,\) which satisfies \(\bar{x}(0)=x_0, \bar{y}(0)=y_0\) and \(\bar{u}(0)=u_0\). According to Remark 1 and 2, this solution satisfies \((\bar{x}(t),\bar{y}(t),\bar{u}(t))\in {\mathbb {R}}^4_+\times \mathcal {U}\) for all \(t\ge 0\).

From Eq. (12), this solution satisfies

$$\begin{aligned} \varphi (\bar{y}(t),h(\bar{x}(t),\bar{y}(t),\bar{u}(t)))\le 0~\hbox {for all}~t\ge 0, \end{aligned}$$

and out of Eq. (7b), it follows that

$$\begin{aligned} \varphi (\bar{y}(t),\dot{\bar{y}}(t))\le 0~\hbox {for all}~t\ge 0. \end{aligned}$$

In addition, since function \(\varphi\) is an \({\mathbb {R}}^m_+\)-Lyapunov function in the sense of definition 2, then we get \(\bar{y}(t)\rightarrow 0\) when \(t\rightarrow \infty\). Consequently, \(\zeta\) provides a feedback control strategy for each \((x_0,y_0,u_0)\in \mathcal {D}_\varphi \cap {\mathbb {R}}^4_+\times \mathcal {U}\). \(\square\)

Now, let us introduce the following set-valued map, for each \(\mu >0,\)

$$\begin{aligned} \mathcal {C}_\varphi ^\mu (x,y,u)\doteq \left\{ v\in \mathcal {G}(x,y,u)~\mid ~ \right. \nonumber \\ \left. \ell _\varphi (x,y,u)+\left\langle m_\varphi (x,y,u),v\right\rangle \le -\mu \right\} , \end{aligned}$$
(21)

for all \((x,y,u)\in \mathcal {K},\) where \(\ell _\varphi\) and \(m_\varphi\) are, respectively, given by (18) and (19). Arguing as in the proof of (Boujallal et al. 2021b, Theorem 2), we can easily get the following result.

Theorem 2

Assume that, for some \(\mu >0\), the map \(\mathcal {C}_\varphi ^\mu\) admits a continuous selection. Then, such selection provides feedback control strategies for each \((x_0,y_0,u_0)\in \mathcal {K}\cap {\mathbb {R}}^4_+\times \mathcal {U}\).

4 Characterization of Continuous Selection

Now, we will characterize the continuous selection of the feedback map \(\mathcal {C}_\varphi ^\mu\), given by (21), and which bring us the vaccination and treatment strategies in order to solve the problem (). For that purpose, let us rewrite the model (3) in the general form of system 7. Let \(x=(x_1,x_2,x_3)^\prime =(S,E,R)^\prime\), \(y=I\) and \(u=(u_1,u_2)^\prime\), then the dynamic

$$\begin{aligned} f(x,y,u)=\begin{pmatrix} \Lambda -c\,\dfrac{x_1y}{N}-dx_1-u_1x_1 \\ c\,\dfrac{x_1y}{N}-(e+d)x_2\\ gy-dx_3+u_1x_1+u_2y \end{pmatrix}, \end{aligned}$$
(22)

and the infectious dynamics is given by

$$\begin{aligned} h(x,y,u)=ex_2-(g+a+d+u_2)y. \end{aligned}$$
(23)

Now, let us express the Lyapunov function \(\varphi\), given in Definition 2, by

$$\begin{aligned} \varphi (y,h(x,y,u))=h(x,y,u)+\alpha y, ~~\hbox {for all}~~\alpha >0. \end{aligned}$$

Indeed, we get thanks to (7b),

$$\begin{aligned} \dot{\bar{y}}(t)\le -\alpha \bar{y}(t), \end{aligned}$$

which yields the estimate

$$\begin{aligned} 0\le \bar{y}(t)\le y_0\exp (-\alpha t) ~~\hbox {for all}~t\ge 0. \end{aligned}$$

Letting \(t\rightarrow \infty\), we get

$$\begin{aligned} \bar{y}(t)\rightarrow 0. \end{aligned}$$

In this case, the expressions of functions \(\ell _\varphi\) and \(m_\varphi\) are given by

$$\begin{aligned} \ell _\varphi (x,y,u)&=ef_2(x,y,u)-(g+a+d+u_2-\alpha )\times \nonumber \\&~~~~h(x,y,u)+\delta u_2y,&\nonumber \\ m_\varphi (x,y,u)&=(0,-y)^\prime .&\end{aligned}$$
(24)

4.1 The Case of Limited Vaccination Constraint

The set of mixed constraints is given as follows:

$$\begin{aligned} \mathcal {K}=\left\{ (x,y,u)\in {\mathbb {R}}^5~\mid ~u_1x_1\le V_0 \right\} . \end{aligned}$$

Let \(\mathcal {D}_\varphi \doteq \left\{ (x,y,u)\in \mathcal {K} ~\mid ~ \phi (x,y,u)\le 0 \right\} ,\) then the maps \(\mathcal {G}\) and \(\mathcal {C}_{\mu }\) are expressed as follows:

$${\mathcal{G}}_{\varphi } (x,y,u) = ~\left| {\begin{array}{*{20}l} {{\mathcal{V}}~{\text{ }}} & {{\text{if }}u_{1} x_{1} < V_{0} ,} \\ {\left\{ {v \in {\mathcal{V}}~{\mid }~u_{1} f_{1} (x,y,u) + x_{1} v_{1} \le 0} \right\}} & {~{\text{ if }}u_{1} x_{1} = V_{0} ,} \\ \end{array} } \right.$$

and

$$\begin{aligned} \mathcal {C}_\varphi ^{\mu }(x,y,u)=\left\{ v\in \mathcal {G}(x,y,u)~\mid ~\ell _\varphi -yv_2\le -\mu \right\} . \end{aligned}$$

Hence, the continuous selection of the map \(\mathcal {C}_{\mu }\) can be expressed as follows:

$$\begin{aligned} \zeta (x,y,u)=(v_1,v_2)^\prime , \end{aligned}$$
(25)

where

$$\begin{aligned} \begin{array}{lll} v_1=\max \left( 0,\min \left( \delta u_1^{\max }, -\dfrac{u_1}{x_1}f_1(x,y,u)\right) \right) ,\\ v_2=\min \left( \delta u_2^{\max },\max \left( 0, \dfrac{\ell (x,y,u)+\mu }{y}\right) \right) . \end{array} \end{aligned}$$

4.2 The Case of Unlimited Vaccination Constraint

The set of constraint is given by

$$\begin{aligned} \mathcal {K}={\mathbb {R}}^5. \end{aligned}$$

Then, the maps \(\mathcal {G}\) and \(\mathcal {C}_\varphi ^{\mu }\) are expressed as follows:

$$\begin{aligned} \displaystyle { \mathcal {G}(x,y,u) = \mathcal {V}}, \end{aligned}$$

and

$$\begin{aligned} \mathcal {C}_\varphi ^{\mu }(x,y,u)=\left\{ v\in \mathcal {V}~\mid ~\ell _\varphi (x,y,u)-yv_2\le -\mu \right\} . \end{aligned}$$

Hence, the continuous selection of the map \(\mathcal {C}_\varphi ^{\mu }\) can be expressed as follows:

$$\begin{aligned} \zeta (x,y,u)=(v_1,v_2)^\prime , \end{aligned}$$
(26)

where

$$\begin{aligned} \begin{array}{lll} &v_1\in [0,\delta u_1^{\max }],\\& v_2=\min \left( \delta u_2^{\max },\max \left( 0, \dfrac{\ell (x,y,u)+\mu }{y}\right) \right) . \end{array} \end{aligned}$$
Fig. 1
figure 1

Simulation results of the model (3) with unlimited vaccination resources (i.e. without mixed constraints). (a) Evolution in time of susceptible population (S) with \(S_0=2\times 10^6\). (b) Evolution in time of exposed population (E) with \(E_0=0.5\times 10^3\). (c) Evolution of infectious population (I) with \(I_0=10^{3}\). (d) Evolution in time of vaccinated population (W) with \(W_0=0\)

Full size image
Fig. 2
figure 2

(a) Evolution of the control \(u_1\) with \(u_1^0=0.4\). (b) Evolution of the control \(u_2\) with \(u_2^0=0\). The control \(u=(u_1,u_2)^\prime\) is given via selection \(\zeta\), given by (26), with \(\mu =2\) and \(\delta =0.8\)

Full size image
Fig. 3
figure 3

Simulation results of the model (3) with limited vaccination resources (i.e. with mixed constraint: \(u_1S\le V_0\) and \(V_0= 10^4\)). (a) Evolution in time of susceptible population (S) with \(S_0=2\times 10^6\). (b) Evolution in time of exposed population (E) with \(E_0=0.5\times 10^3\). (c) Evolution of infectious population (I) with \(I_0=10^{3}\). (d) Evolution in time of vaccinated population (W) with \(W_0=0\)

Full size image
Fig. 4
figure 4

(a) Evolution of the control \(u_1\) with \(u_1^0=5\times 10^{-3}\). (b) Evolution of the control \(u_2\) with \(u_2^0=0\). The control \(u=(u_1,u_2)^\prime\) is given via selection \(\zeta\), given by (25), with \(\mu =2\) and \(\delta =0.8\)

Full size image
Fig. 5
figure 5

A comparison of the evolution in time of susceptible, exposed and infected population without and with mixed constraint (with \(V_0= 10^4\)).(a) Evolution in time of susceptible population (S) with \(S_0=2\times 10^6\). (b) Evolution in time of exposed population (E) with \(E_0=0.5\times 10^3\). (c) Evolution of infectious population (I) with \(I_0=10^{3}\)

Full size image
Fig. 6
figure 6

Simulation results of the model (3) with limited vaccination resources and different values of \(V_0\). (a) Evolution in time of susceptible population (S). (b) Evolution in time of exposed population (E). (c) Evolution of infectious population (I). (d) Evolution in time of vaccinated population (W)

Full size image

5 Numerical Simulation

In this section, numerical simulations of the system (3) are carried out in order to describe the spread of COVID-19 with treatment and under the limitation of vaccination resources. All simulations are performed using MATLAB and the parameters are given in Table 1. It should be noted that inspired by Carcione et al. (2020), the dynamic of the vaccinated individuals, denoted by W, is given by \(\dot{W}=u_1S\).

Table 1 Parameter’s description and value for model (1)
Full size table

In Fig. 1, we consider the case where the two \(u_1\) and \(u_2\) controls are both applied and where there are no mixed constraints. The Fig. 1a–d show the number of people in the compartments (S), (E), (I) and (W) with and without controls. Our numerical results show that the use of both controls significantly reduces the number of susceptible, exposed and infected people. This significant decline can be explained by the availability of vaccination resources (the absence of mixed vaccination constraint) which allows to vaccinate a very large number of the population and save the lives of more than 786706 people (see Fig. 1d). In this case, the controls \(u_1\) and \(u_2\) are given in Fig. 2.

For Fig. 3, it concerns the case where the supply of vaccines is limited with \(u_1S \le V_0 = 10^4\). It can be seen in Fig. 3b, c that even with constraints on the number of people vaccinated; vaccination and treatment combined can help reduce the number of people exposed and infected with COVID-19. On the other hand, the limit of resources available for a good vaccination campaign means that the number of susceptible or vulnerable people remains higher compared to the case without control (see Fig. 3a). Also, we note in Fig. 3d that the number of people vaccinated will not exceed 12,500 people, i.e. a drop of 98.41% compared to the result obtained in Fig. 1d. Profiles of the controls \(u_1\) and \(u_2\) are depicted in Fig. 4.

Figure 5 shows the evolution of the number of individuals in compartments (S), (E) and (I) in both cases: with and without mixed constraints. It can be seen that the numbers of individuals in these compartments remain higher in the case where mixed constraints are considered.

In Fig. 6, we depict the dynamic of our controlled model using different values of \(V_0\). In this figure, we can see the important impact of \(V_0\) on the evolution of the number of susceptible, exposed, infected and vaccinated individuals. More \(V_0\) is important more the number of people in S, E and I compartments decrease. So, it is important to double the global efforts to give access to vaccines, especially for low-income countries, which can help protect a maximum number of people against the SARS-CoV-2 virus.

6 Concluding Remarks and Future Works

Through this work, we contribute to the global efforts that aim to understand and control the spread of COVID-19 disease. We considered a compartmental epidemic model where the population is divided into four categories (susceptible, exposed, infectious and removed individuals) and includes two control variables, namely vaccination and treatment. As it is well known, vaccination against COVID-19 provides multiple challenges that are not encountered by many other vaccination programmes. In several countries, infection rates are increasing exponentially, while vaccination rates are restricted by supply and logistics. In this context, we formulated a control problem with mixed constraints that takes into account the limitation of vaccination resources and aims to reduce asymptotically the number of infected people to zero. To solve our control problem, we applied a set-valued approach based on viability theory and set-valued analysis. More precisely, the control pair that reduces the number of infectious (I) individuals is introduced via continuous selections of a designed feedback map.

In future work, it will be interesting to extend this study to other diseases and consider other mixed constraints related to age and regional spread of the disease.