A direct integral pseudospectral method for solving a class of infinite-horizon optimal control problems using Gegenbauer polynomials and certain parametric maps

1.   Introduction

Nowadays, mathematical modeling has become an essential tool in the fight against infectious diseases and epidemics. A good model serves to describe the evolution of the epidemic, to provide possible solutions, to predict the asymptotic behavior, to succeed in computer simulations, etc. ^[1,2,3,4,5].

In the field of disease modeling, numerous approaches have been developed to understand and predict disease evolution. These models play a crucial role in understanding the dynamics of infectious diseases and in guiding decision-making processes (see ^{[6,7,8,9,10,11,12,13]}). By capturing the complex interplay between factors such as population dynamics, transmission mechanisms and intervention strategies, these models enable researchers to gain a comprehensive understanding of disease progression.

So scientists everywhere in the world, are always looking for a perfect mathematical model. Our contribution to this theme consists of adopting the $SIR$ compartment model described by the following nonlinear differential equations:

where $S, \, I\:$and $R\,$ are respectively the densities of susceptible individuals, infectious individuals and recovered individuals, $[0, T_{f}]\,$ is the duration of time over which we study our system, $\varLambda\,$ is the Susceptible individual recruitement, $m$ is the natural mortality rate, $r$ is the mortality rate due to the virus, $g$ is the natural recovery rate and $\alpha\,$ is the disease transmision coefficient.

Of course, the coefficient $\alpha$, which indicates the frequency at which susceptible people contract the infection, is the most important parameter of the model; therefore, faulty knowledge of the latter leads to a model that does not reflect the evolution of the disease at all. To justify this we give as an example the following simulations to show that a slight modification of the transmission rate $\alpha$ induces a large difference in the number of infected individuals. Indeed, in the following figure, we establish the curves representing the number of infected individuals for $\alpha_{1} = 0.3, \alpha_{2} = 0.4$ and $\alpha_{3} = 0.5$.

As indicated in the abstract of this paper, our main objective is to identify the parameter $\alpha.$ Indeed, the third and the fourth sections of our work will be devoted to the following inverse problem : Knowing the number of infected individuals I(t), we are inspired by the works of ^[14,16] to reconstruct the parameter $\alpha$.

The strategy we will adopt is as follows, for $0 < T^{\ast} < T_{f}\,$we consider the model

and, we use the output equation $f(t) = I(t)$ (prevalence function) to determine the transmission rate function $\alpha(t)$. Once $\alpha(t)$ is found, we take the arithmetic mean $\alpha_{m} = \frac{1}{T^{\ast}}\int_{0}^{T^{\ast}}\alpha(t)dt$ (this mean is calculated by using MATLAB). In the fourth section, and to take advantage of the previous results of this paper, we start from the valid module (with the identified $\alpha$) and we are interested in the following control problem

with

where the controls $u$ and $v$ mean respectively the vaccine and the treatment applied to the susceptible and infected persons. It is clear that the optimization of $J$ implies that we are interested in the controls that achieve the following objectives:

1) The number of infected individuals $I(t)\,$ must be minimal during the period $[T^{\ast}, T_{f}]$.

2) The number of recoveries $R(t)\,$ should be as high as possible during the period $[T^{\ast}, T_{f}]$.

3) The costs $\parallel u^{\ast}\parallel\,$ and $\parallel v^{\ast}\parallel$ must be very minimal.

4) $T^{\ast}\,$ must be minimal, which means that the duration of identification must be as short as possible.

It should be noticed that the fundamental result of this work, namely, the identification of the infection rate $\alpha$ allows us to achieve the following goal:

To have a reliable mathematical model that can be used for decision support, such as determining a good control measure (subject of Section 5), predicting the peaks of the epidemic, determining the equilibrium points and studying their stability, etc.

When South Africa declared to the World Health Oraganization (WHO) the appearance of the new Omicron variant of COVID-19, the international community took many measures that seriously affected the economy of South Africa (air embargo, cessation of a product port from South Africa, etc.) ^[17]. Therefore, this will encourage many countries (especially countries with weak economies) to no longer declare their real health situation to the WHO.

Our work allows us then to bring a solution to such situations, and as we have seen during these four years during the COVID-19 pandemic, every country must declare the number of infected individuals "$I(t)$" to the WHO, to take advantage of the financial and logistic help. Thus, using the function of prevalence of $I(t)$, the (WHO) can detect any new variants of COVID-19 (see the last figure of Section 3).

The structure of this paper is as follows: In Section 2, we prove the existence, the positivity and the boundedness for our system. In Section 3, we demonstrate the identification of the transmission rate and we apply it to real data concerning Morocco. The existence of an optimal solution and the neccesary optimality conditions are confirmed in Section 4. Regarding application, the numerical results related to our control problem are given in Section 5. In the end, we conclude the article in Section 6.

2.   Model basic properties

As indicated in the introduction of this paper, we consider the SIR compartment model defined by

with

All of the system's parameters have the meanings presented in Table 1.

2.1.   Positivity and boundedness of solutions

For biological reasons and inspired by ^{[18,19,20,21]}, we establish the the following result.

Theorem 1. With a positive initial condition, the solution $(S(t), \; I(t), \; R(t))$ of (2.1) is positive and eventually uniformly bounded on $[0, +\infty)$.

Proof. Assume that the solution $(S(t), \; I(t), \; R(t))$ with a positive initial condition exists and is unique on $[0, b)$, where $0 < b < +\infty$

we have

Similarly,

Integrating the above inequality from 0 to $t$ yields

using the same approach, we can demonstrate that $R(t) > 0\, \;$for all $t\in[0, b).$ Therefore, we obtain that $S(t) > 0, \; I(t) > 0$ and $R(t) > 0$ for all $t\in[0, b)$; furthermore, $N(t)$ the total population studied satisfies the following equation

which implies that

thus $(S(t), \; I(t), \; R(t))$ is bounded on $[0, b)$; therefore, we have that $b = +\infty$, $\underset{t\rightarrow +\infty}{lim}supN(t)\leq\frac{\Lambda}{m}.$

Hence $(S(t), \; I(t), \; R(t))$ are all identically bounded by a fixed upper bound $\frac{\Lambda}{m}$. □

2.2.   Existence of the solution

Theorem 2. For given initial conditions $S(0)\geq0, \; I(0)\geq0$ and $R(0)\geq0$ the system (2.1) admits a unique solution.

Proof. Let $\phi(t) = \left(\begin{array}{c} S(t)\\ I(t)\\ R(t) \end{array}\right)$ and $\phi_{t}(t) = \left(\begin{array}{c} \frac{dS}{dt}\\ \frac{dI}{dt}\\ \frac{dR}{dt} \end{array}\right)$;

so the system can be rewritten in the following form

where

$A = \left(\begin{array}{ccc} -m & 0 & 0\\ 0 & -(m+r+g) & 0\\ 0 & g & -m \end{array}\right)\,$ and $N(\phi) = \left(\begin{array}{c} \Lambda-\alpha S(t)I(t)\\ \alpha S(t)I(t)\\ 0 \end{array}\right)$,

and $\phi_{t}$ denotes the derivative of $\phi$ with respect to time.

The second term on the right-hand side satisfies

by putting $M = max(2k\underset{t\in[0, T]}{sup}\alpha(t)) < \infty$

thus it follows that the function $N$ is uniformly Lipshitz-continuous; therefore, the solution of the system exists. □

3.   Identification of the transmission rate

In this section, we start with a database containing the number of infected people daily $I(t_{1}), I(t_{2}), ...., I(t_{n})$, where "$t_{i}$" designates the day "$i$" (this database, which contains data pertaining to each country, was produced by JavaScript Object Notation (JSON) ^[22]. Then, we use a polynomial interpolation to obtain a function "$I(t)$", and the interpolation polynomial $I(t)$ can be used to estimate the number of infected persons at each time $t\,$ belonging to $[0, T]$, where $T$ is the duration of the epidemic in the country under consideration (for COVID-19, $T$ = 4 years and $I(t)\,$is obtained by using Python).

This function $I(t)$ reffered to as the prevalence function, is denoted as $f(t)$ is at the origin of the following fundamental result.

Theorem 3. Given $m, \; g, \; r > 0, \; \alpha_{0} > 0$ and $T > 0,$ there exists $K > 0$ such that if $\alpha_{0} < K$ there is a solution $\alpha(t)$ with $\alpha(0) = \alpha_{0}$ such that $I_{\alpha}(t) = f(t)$ for $0\leq t\leq T$, where $\left(\begin{array}{c} S_{\alpha}(.)\\ I_{\alpha}(.)\\ R_{\alpha}(.) \end{array}\right)$ is the solution of system (2.1) corresponding to $\alpha(t)$ if $\frac{f'(t)}{f(t)} > -(m+r+g)$ for $0\leq t\leq T.$

Proof. The second equation of system (2.1) implies that,

which must be positive for $0\leq t\leq T$.

By rewriting (3.1) as

we can compute $S'(t)$ and substitute it for the first equation of system (2.1) in order to find that

which will provide for us later

We set

and

Thus, our equation is

which gives us the following

The nonlinear ODE can be converted to linear by setting $x(t) = \frac{1}{\alpha(t)},$ which means that $x'(t) = \frac{-\alpha'(t)}{\alpha(t)^{2}}$

We know that the explicit solution is provided by the approach of integrating factors

where $z(t) = \int_{0}^{t}z(s)ds$.

Then, $\alpha(0)$ needs to satisfy certain conditions in order for $\alpha(t)$ to be positive:

So according to the mathematics, there is an infinite number of possible values for $\alpha(0)$, and as a result, an infinite number of transmission functions $\alpha(t)$.

The inverse problem is thus inadequately characterized

□

3.1.   Application (case of Morroco)

We will present a mathematical approach showing how to detect a change in the veracity of the virus (i.e., the appearance of a possible mutation) in any country based on official data. In our work, we will focus on the daily infection data of COVID-19 in Morocco.

Below is a representative figure of the daily infection of COVID-19 in Morocco since the appearance of the virus in early March 2020 and up to 700 days later ^[23].

In order to extract the transmission rate of COVID-19 in Morocco from the simulated data, we will use a polynomial interpolation. More precisely, we consider the data of 700 days from the first day of appearance of COVID-19 in Morocco and which appears in the previous table (Figure 2); we divide the data of 700 days into five periods $P_{1}, P_{2}, P_{3}, P_{4}, P_{5}$ each comprising 140 days. Then, for the daily infections recorded for the 140 days of each period $P_{i}$, we take the arithmetic mean $I_{i}$. Thus and starting from the table of the averages $I_{1}, I_{2}, I_{3}, I_{4}, I_{5},$ we construct the Lagrange interpolation polynomial $P(t)$. In order to succeed in this task and in view of the complexity of the calculations, we chose to use Python code to establish the interpolation polynomial $P(t).$ To activate Python, we developed the following appropriate code:

Therefore, using Python, we can construct the interpolation polynomial

Thus, we have constructed a function $I(t) = P(t)$ which allows us to estimate the number of infected individuas at any time point within the period $[March\, 3, 2020, March\, 3, 2020\, +700\, day]$.

Consequently, we inject the function $I(t)$ given by (3.9) ($I(t) = P(t)$) in Theorem 3 to obtain the function of transmission $\alpha(t)$ whose curve is generated via MATLAB, (see Figure 3).

Through the two previous figures, we show the validity of the results of Theorem 3, and particularly for the case of Morocco. Indeed, from these two figures, we notice that over the same period of time $[t_{100}, t_{220}]$, the number of infected people $I(t)\,$ also increased as the value of $\alpha(t)$ increased, and this period corresponds perfectly to the appearance of the Delta variant in Morroco. Similarly, over the period $[t_{550}, t_{700}]$, another increase of infected people occurred as the transmission function increased, and this is exactly the date of the appearance of the Omicron variant in Morocco.

Therefore, to detect new variants of COVID-19 or any virus it is sufficient to monitor and analyze the infection rate of individuals within a region or country.

Remark 4. Having the database corresponding to the daily data of the infected in a given country for a period of $N$ days, (or during the time horizon $[0, T]$), if we are interested in identifying $\alpha(t)$ only over a period of time $[0, T^{\ast}]\subseteq[0, T]$ (i.e., $T^{\ast} < T$), we can easily adapt the previous code. To better explain the objective of this remark, we take the case of Morocco and expose the modification of the previous PYTHON code to $\alpha(t)$ only based on the interpolation polynomial corresponding to a period $T^{\ast} = 300$ and $n = 4$ ($n$ is the chosen degree of the interpolator polynomial).

4.   An optimal control problem

In this section, we are interested in the $SIR$ uncontrolled model given by

where $S, \; I$ and $R$ are the densities of susceptible, infected and removed people respectively, and the transmission coefficient $\alpha$ is assumed to be unknown. We assume to have access to $I(t)$ (i.e the number of infected individuals) at each time $t$. Our objective is to determine the control laws $u(t)$ and $v(t)$ ($u(t)$ is the vaccine applied to susceptible persons $S(t)$ and $v(t)$ is the treatment applied to infected persons $I(t)$) that will allow us to reach the following objectives:

and do this with minimal effort and cost. Mathematically, this can be interpreted by the optimization of the function

under the following constraints

But as mentioned at the beginning of this section the transmission coefficient $\alpha$ is not known; therefore, the model (4.3) is not reliable. To face this problem, we will segment our time interval in the following way.

Here the durability $[0, T^{\ast}]$ is devoted to the identifiability of the unknown parameter $\alpha$ ($T^{\ast}$ is as low as possible), and the horizon time $[T^{\ast}, T_{f}]$ is devoted to the realization of previous objectives. Indeed, to do that we adapt the following strategy

First, we use the prevalence function $I(t), \; t\in[0, T^{\ast}]$, which has been obtained by interpolation of the database corresponding to $T^{\ast}$ (see Remark 4, Section 3).

Second, we use Theorem 3 to identify the transmission function $\alpha(t)$ corresponding to the following model:

Third, we assign to our unknown parameter $\alpha$ the mean value $\alpha_{m} = \frac{1}{T^{\ast}}\int_{0}^{T^{\ast}}\alpha(t)dt$, which is calculated by using MATLAB (for the Morrocan case, $T^{\ast} = 60\; days$ the Python code gives the prevalence function $I(t)$ which is injected in Theorem 3 to obtain via Matlab the transmission function $\alpha(t)$ for the horizon time $[0, T^{\ast}]$ and consequently $\alpha_{m} = 0.011$).

Remark 5. We chose $T^{\ast} = 60$, which led us to $\alpha_{m} = 0, 011$ and which coincides exactly with the transmission coefficient used in the COVID-19 model adopted by the Moroccan authorities ^[24], because we noticed that from the beginning of the third month of the epidemic the numbers of infected and recovered people stabilized. The following table shows that the value of $\alpha_{m} = 0.011$ is indeed admissible.

4.1.   Statement of the problem

Having the value of $\alpha_{m} = 0.011$, we use the well identified model,

with the initial given data given by

to reach our predefined objectives ^[25], i.e., to minimize the cost function:

where $\rho_{i}, \; i = 1, 2, 3$ are are weightings of the positive factor parameters associated with the controls $u$ and $v$ and the time, respectively. The presence of the term $\rho_{3}(T_{f})^{2}$ in (4.3) of the cost function means that we introduce another objective in our project, namely: The search for the controls $u$ and $v$ that allow one to minimize the number of infected people, maximize the number of recovered, and to do this with a minimal cost and during a period $T_{f}$ as short as possible.

In other words, we search for the optimal controls $u^{\ast}$ and $v^{\ast}$ and the optimal time $T_{f}$ such that

Subject to system (4.5), we take $U_{ad}$ as the control set defined by

4.2.   The existence of the optimal solution

In this section, we use the classical results of Fleming (1975) and Luckes (1982) ^[26,27] to prove the existence of an optimal control measure. Effectively, let $f = [f_{1}, f_{2}, f_{3}]^{\intercal}$ be the right-hand side of the system (4.5) and $X = [S, I, R]^{\intercal}$.

Theorem 6. There exists control functions $u^{\ast}$ and $v^{\ast}$ so that

Proof. It is simple to confirm that an optimal control exists in order to demonstrate the following:

(A1) The set of controls and corresponding state variables is not empty.

(A2) The admissible set $U_{ad}$ is convex and closed.

(A3) A linear function comprising the state and control variables bounds the right side of the state system.

(A4) The integrand of the objective function is convex on $U_{ad}.$

(A5) There exists constants $\omega_{1} > 0, \; \omega_{2} > 0$ and $\psi > 1$ such that the integrand $L(I, R, u, v)$ of the objective function satisfies

(A6) The pay-off term at the optimal time in the objective function $\varphi(X(T_{f}), T_{f}) = \rho_{3}(T_{f})^{2}.$

We use a simplified version of an existence result to prove that the set of controls and corresponding state variable is not empty for the first condition ^[28] (see Theorem 7.1.1).

The set $U_{ad}$ is convex and closed by definition; thus, the condition (A2) holds.

Our state system is bilinear in $u$ and $v$. Moreover, the solutions of the system are bounded; hence, the condition (A3) holds.

Note that the integrand of our objective function is convex from where condition (A4) is verified.

To prove the condition (A5), we have that $u+v\geq\mid u\mid^{2}+\mid v\mid^{2}$ since $0\leq u, v\leq1$.

It is obvious that $\varphi(X(T_{f}), T_{f})$ is continuous; hence, the last condition is required.

The result follows directly from ^[26]. □

4.3.   Characterization of the optimal control scheme

It is easy to establish that the corresponding Lagrangian $L$ and Hamiltonian $H$ are, respectively,

with

and where $\lambda_{i}, i = 1, 2, 3$ are the adjoint functions to be determined suitably.

Next, by applying the Pontryagin maximum principle ^[29] to the Hamiltonian $H$, we get the following theorem.

Theorem 7. Given the optimal control scheme $(u^{\ast}, v^{\ast})$, the optimal time $T^{\ast}$ and the solutions $S^{\ast}, \; I^{\ast}\; and\; R^{\ast}$ of the corresponding state system (4.5), there exists adjoint variables $\lambda_{1}, \; \lambda_{2}$ and $\lambda_{3}$ satisfying

with the transversality condition $\lambda_{i}(T_{f}) = 0$ for $i = 1, 2, 3.$

Furthermore, the optimal control scheme $(u^{\ast}, v^{\ast})$ is given by

and the optimal time is given by

Proof. By applying Pontryagin's maximum principle to the state, we obtain the adjoint equations of the problem by direct computation

and the controls are obtained by calculating the following equations

with the transversality conditions

By the bounds in $U_{ad}$ of the controls ^[30], it is easy to obtain that $u^{\ast}$ and $v^{\ast}$ are given in the form (4.17) of system (4.5).

The transeversality condition for $T_{f}^{\ast}$ to be the optimal time can be stated as

where

then,

□

5.   Numerical simulations

In this section, we present numerical results that clearly illustrate the reliability of our approach. For this purpose, we have opted for the Runge Kutta method of order 4 and we have used the initial data conditions $T^{\ast} = 60$, $S(T^{\ast})$, $I(T^{\ast})$ and $R(T^{\ast})$ from our database and a MATLAB program to obtain the following results.

In Table 3, we present the values of the initial conditions, parameters and constants, which have been used in our computations ^[31].

Figures 4 and 5 show that our obtained optimal controls allows us to clearly reduce the number of susceptibles and infected people.

In Figure 6, we see that the final optimal time corresponding to our problem is $T_{f} = 245$ days, which was derived from the curve function $G(t) = H(t)+\frac{\partial\varphi(t)}{\partial t}$.

Consequently, the advantages of our method stem from the combined use of real data, customized modeling and optimal control. Based on country-specific data, our simulations capture the unique characteristics of the local epidemic, enabling more accurate prediction and control strategies. The incorporation of optimal control techniques enables us to optimize interventions, taking into account various constraints and objectives.

6.   Conclusions

In this work, we were interested in two fundamental aspects of systems theory, namely identification and regulation. Indeed, we considered a classical SIR model corresponding to the evolution of COVID-19 in a given country and where the transmission coefficient is supposed to be unknown. We were inspired by the work on inverse problems developed in ^[15,16] to determine the infection coefficient, using the starting point of the prevalence function $I(t)$ = number of infected at time $t$. We made use of a country-specific database and the features of the programming languages MATLAB and Python to complete this work successfully. After identifying the mathematical model, we used an optimal control problem to address the difficulties given by the identification procedure. As an illustration, our study explicitly looked at the instance of Morocco. We made great progress in solving several COVID-19-related issues with this effort. To further increase the mathematical model's applicability and effectiveness, it is imperative to emphasize the importance of understanding all four parameters. We can offer more thorough answers and insights into how to handle the complexity of COVID-19 by accepting these developments;

This will be the subject of future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Kamal Shah and Thabet Abdeljawad are thankful to Prince Sultan University for paying the APC and for the support through the TAS research lab.

Conflict of interest

The authors declare that they have no conficts of interest.