
In this work, we formulate a fractal fractional chaotic system with cubic and quadratic nonlinearities. A fractal fractional chaotic Lorenz type and financial systems are studied using the Caputo Fabrizo (CF) fractal fractional derivative. This study focuses on the characterization of the chaotic nature, and the effects of the fractal fractional-order derivative in the CF sense on the evolution and behavior of each proposed systems. The stability of the equilibrium points for the both systems are investigated using the Routh-Hurwitz criterion. The numerical scheme, which includes the discretization of the CF fractal-fractional derivative, is used to depict the phase portraits of the fractal fractional chaotic Lorenz system and the fractal fractional-order financial system. The simulation results presented in both cases include the two- and three-dimensional phase portraits to evaluate the applications of the proposed operators.
Citation: Ihtisham Ul Haq, Shabir Ahmad, Sayed Saifullah, Kamsing Nonlaopon, Ali Akgül. Analysis of fractal fractional Lorenz type and financial chaotic systems with exponential decay kernels[J]. AIMS Mathematics, 2022, 7(10): 18809-18823. doi: 10.3934/math.20221035
[1] | Asma Hanif, Azhar Iqbal Kashif Butt, Tariq Ismaeel . Fractional optimal control analysis of Covid-19 and dengue fever co-infection model with Atangana-Baleanu derivative. AIMS Mathematics, 2024, 9(3): 5171-5203. doi: 10.3934/math.2024251 |
[2] | Aqeel Ahmad, Cicik Alfiniyah, Ali Akgül, Aeshah A. Raezah . Analysis of COVID-19 outbreak in Democratic Republic of the Congo using fractional operators. AIMS Mathematics, 2023, 8(11): 25654-25687. doi: 10.3934/math.20231309 |
[3] | Moh. Mashum Mujur Ihsanjaya, Nanang Susyanto . A mathematical model for policy of vaccinating recovered people in controlling the spread of COVID-19 outbreak. AIMS Mathematics, 2023, 8(6): 14508-14521. doi: 10.3934/math.2023741 |
[4] | Sohail Ahmad, Ponam Basharat, Saleem Abdullah, Thongchai Botmart, Anuwat Jirawattanapanit . MABAC under non-linear diophantine fuzzy numbers: A new approach for emergency decision support systems. AIMS Mathematics, 2022, 7(10): 17699-17736. doi: 10.3934/math.2022975 |
[5] | Shuyun Jiao, Mingzhan Huang . An SIHR epidemic model of the COVID-19 with general population-size dependent contact rate. AIMS Mathematics, 2020, 5(6): 6714-6725. doi: 10.3934/math.2020431 |
[6] | Deniz UÇAR, Elçin ÇELİK . Analysis of Covid 19 disease with SIR model and Taylor matrix method. AIMS Mathematics, 2022, 7(6): 11188-11200. doi: 10.3934/math.2022626 |
[7] | Yilihamujiang Yimamu, Zui-Cha Deng, Liu Yang . An inverse volatility problem in a degenerate parabolic equation in a bounded domain. AIMS Mathematics, 2022, 7(10): 19237-19266. doi: 10.3934/math.20221056 |
[8] | Khalaf M. Alanazi . The asymptotic spreading speeds of COVID-19 with the effect of delay and quarantine. AIMS Mathematics, 2024, 9(7): 19397-19413. doi: 10.3934/math.2024945 |
[9] | CW Chukwu, S. Y. Tchoumi, Z. Chazuka, M. L. Juga, G. Obaido . Assessing the impact of human behavior towards preventative measures on COVID-19 dynamics for Gauteng, South Africa: a simulation and forecasting approach. AIMS Mathematics, 2024, 9(5): 10511-10535. doi: 10.3934/math.2024514 |
[10] | Qinyue Zheng, Xinwei Wang, Qiuwei Pan, Lei Wang . Optimal strategy for a dose-escalation vaccination against COVID-19 in refugee camps. AIMS Mathematics, 2022, 7(5): 9288-9310. doi: 10.3934/math.2022515 |
In this work, we formulate a fractal fractional chaotic system with cubic and quadratic nonlinearities. A fractal fractional chaotic Lorenz type and financial systems are studied using the Caputo Fabrizo (CF) fractal fractional derivative. This study focuses on the characterization of the chaotic nature, and the effects of the fractal fractional-order derivative in the CF sense on the evolution and behavior of each proposed systems. The stability of the equilibrium points for the both systems are investigated using the Routh-Hurwitz criterion. The numerical scheme, which includes the discretization of the CF fractal-fractional derivative, is used to depict the phase portraits of the fractal fractional chaotic Lorenz system and the fractal fractional-order financial system. The simulation results presented in both cases include the two- and three-dimensional phase portraits to evaluate the applications of the proposed operators.
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:
{dSdt=Λ−αS(t)I(t)−mS(t),dIdt=αS(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t),t∈[0,Tf], |
where S,Iand R are respectively the densities of susceptible individuals, infectious individuals and recovered individuals, [0,Tf] is the duration of time over which we study our system, Λ 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 α is the disease transmision coefficient.
Of course, the coefficient α, 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 α 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 α1=0.3,α2=0.4 and α3=0.5.
As indicated in the abstract of this paper, our main objective is to identify the parameter α. 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 α.
The strategy we will adopt is as follows, for 0<T∗<Tfwe consider the model
{dSdt=Λ−α(t)S(t)I(t)−mS(t),dIdt=α(t)S(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t),t∈[0,Tf], |
and, we use the output equation f(t)=I(t) (prevalence function) to determine the transmission rate function α(t). Once α(t) is found, we take the arithmetic mean αm=1T∗∫T∗0α(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 α) and we are interested in the following control problem
J(u∗,v∗,Tf)=min{J(u(t),v(t))/(u,v)∈UadandT∈R∔}, |
with
{dSdt=Λ−αmS(t)I(t)−mS(t)−u(t)S(t),dIdt=αmS(t)I(t)−mI(t)−rI(t)−gI(t)−v(t)I(t),dRdt=gI(t)−mR(t)+u(t)S(t)+v(t)I(t),t∈[T∗,Tf], |
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∗,Tf].
2) The number of recoveries R(t) should be as high as possible during the period [T∗,Tf].
3) The costs ∥u∗∥ and ∥v∗∥ must be very minimal.
4) T∗ 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 α 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.
As indicated in the introduction of this paper, we consider the SIR compartment model defined by
{dSdt=Λ−α(t)S(t)I(t)−mS(t),dIdt=α(t)S(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t),t∈[0,Tf], | (2.1) |
with
S(0)≥0,I(0)≥0andR(0)≥0. | (2.2) |
All of the system's parameters have the meanings presented in Table 1.
Parameter | Signification |
Λ | Susceptible person recruitement |
α(t) | Disease transmission coefficient |
g | Natural recovery rate |
m | Natural mortality rate |
r | Mortality rate due to the virus |
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,+∞).
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<+∞
dSdt=Λ−α(t)S(t)I(t)−mS(t)≥−α(t)S(t)I(t)−mS(t); |
we have
S(t)≥S(0)exp(−∫t0(α(θ)I(θ)+m)dθ)>0. |
Similarly,
dIdt=α(t)S(t)I(t)−(m+g+r)I(t)>−(m+g+r)I(t). |
Integrating the above inequality from 0 to t yields
I(t)≥I(0)exp(−∫t0(m+g+r)dθ)>0; |
using the same approach, we can demonstrate that R(t)>0for all t∈[0,b). Therefore, we obtain that S(t)>0,I(t)>0 and R(t)>0 for all t∈[0,b); furthermore, N(t) the total population studied satisfies the following equation
dNdt=Λ−mN(t)−rI(t)<Λ−mN(t), |
which implies that
N(t)≤Λm+N(0)exp(−mt), | (2.3) |
thus (S(t),I(t),R(t)) is bounded on [0,b); therefore, we have that b=+∞, limt→+∞supN(t)≤Λm.
Hence (S(t),I(t),R(t)) are all identically bounded by a fixed upper bound Λm.
Theorem 2. For given initial conditions S(0)≥0,I(0)≥0 and R(0)≥0 the system (2.1) admits a unique solution.
Proof. Let ϕ(t)=(S(t)I(t)R(t)) and ϕt(t)=(dSdtdIdtdRdt);
so the system can be rewritten in the following form
ϕt(t)=Aϕ+N(ϕ), |
where
A=(−m000−(m+r+g)00g−m) and N(ϕ)=(Λ−αS(t)I(t)αS(t)I(t)0),
and ϕt denotes the derivative of ϕ with respect to time.
The second term on the right-hand side satisfies
∣N(ϕ1)−N(ϕ2)∣=∣N1(ϕ1)−N1(ϕ2)∣+∣N2(ϕ1)−N2(ϕ2)∣=∣α(t)(S1I1−S2I2)∣+∣α(t)(S1I1−S2I2)∣≤2kα(t)(∣S1−S2∣+∣I1−I2∣); |
by putting M=max(2ksupt∈[0,T]α(t))<∞
∣N(ϕ1)−N(ϕ2)∣≤M(∣S1−S2∣+∣I1−I2∣)≤M∥ϕ1−ϕ2∥, |
thus it follows that the function N is uniformly Lipshitz-continuous; therefore, the solution of the system exists.
In this section, we start with a database containing the number of infected people daily I(t1),I(t2),....,I(tn), where "ti" 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,α0>0 and T>0, there exists K>0 such that if α0<K there is a solution α(t) with α(0)=α0 such that Iα(t)=f(t) for 0≤t≤T, where (Sα(.)Iα(.)Rα(.)) is the solution of system (2.1) corresponding to α(t) if f′(t)f(t)>−(m+r+g) for 0≤t≤T.
Proof. The second equation of system (2.1) implies that,
f′(t)+(m+r+g)f(t)=α(t)S(t)f(t), | (3.1) |
which must be positive for 0≤t≤T.
By rewriting (3.1) as
S(t)=f′(t)+(m+r+g)f(t)α(t)f(t), | (3.2) |
we can compute S′(t) and substitute it for the first equation of system (2.1) in order to find that
ddt(f′(t)+(m+r+g)f(t)α(t)f(t))=Λ−α(t)(f′(t)+(m+r+g)f(t)α(t)f(t))f(t)−m(f′(t)+(m+r+g)f(t)α(t)f(t)), | (3.3) |
which will provide for us later
[f"(t)f(t)−f′(t)2+mf(t)(f′(t)+(m+r+g)f(t))]α(t)−[f(t)(f′(t)+(m+r+g)f(t))]α′(t)+[f′(t)+(m+r+g)f(t)−Λ]f(t)2α(t)2=0. |
We set
w(t)=f"(t)f(t)−f′(t)2+mf(t)(f′(t)+(m+r+g)f(t))f(t)(f′(t)+(m+r+g)f(t)), |
and
z(t)=(f′(t)+(m+r+g)f(t)−Λ)f(t)2f(t)(f′(t)+(m+r+g)f(t)). |
Thus, our equation is
α′(t)−w(t)α(t)−z(t)α(t)2=0, | (3.4) |
which gives us the following
α′(t)α(t)2−w(t)α(t)−1−z(t)=0. | (3.5) |
The nonlinear ODE can be converted to linear by setting x(t)=1α(t), which means that x′(t)=−α′(t)α(t)2
x′(t)+w(t)x(t)+z(t)=0. | (3.6) |
We know that the explicit solution is provided by the approach of integrating factors
x(t)=1α(t)=x(0)e−w(t)−e−w(t)∫t0ew(s)z(s)ds, | (3.7) |
where z(t)=∫t0z(s)ds.
Then, α(0) needs to satisfy certain conditions in order for α(t) to be positive:
∫t0ew(s)z(s)ds<1α(0). |
So according to the mathematics, there is an infinite number of possible values for α(0), and as a result, an infinite number of transmission functions α(t).
The inverse problem is thus inadequately characterized
α(t)=1/[e−w(t)α(0)−e−w(t)∫t0ew(s)z(s)ds]. | (3.8) |
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 P1,P2,P3,P4,P5 each comprising 140 days. Then, for the daily infections recorded for the 140 days of each period Pi, we take the arithmetic mean Ii. Thus and starting from the table of the averages I1,I2,I3,I4,I5, 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:
Algorithm 1 Algorithm for extracting the Lagrange interpolation polynomial |
import pandas as pd |
import json |
import matplotlib.pyplot as plt |
import matplotlib.dates as mdates |
import numpy as np |
from scipy.interpolate import lagrange |
from numpy.polynomial.polynomial import Polynomial |
# Opening JSON file |
f = open('country_newcases.json') |
# returns JSON object as # a dictionary |
data = json.load(f) |
df = pd.DataFrame(data) |
# print(df) |
fig, ax = plt.subplots(1, 1, figsize = (10, 7), constrained_layout = True) |
# common to all three: |
ax.plot('last_updated', 'new_infections', data = df) |
# Major ticks every half year, minor ticks every month, |
ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth = (1, 8))) |
ax.xaxis.set_minor_locator(mdates.MonthLocator()) |
ax.grid(True) |
ax.set_ylabel(r'Cases number') |
for label in ax.get_xticklabels(which = 'major'): |
label.set(rotation = 30, horizontalalignment = 'right') |
ax.ticklabel_format(style = 'plain', useOffset = False, axis = 'y') |
ax.set_title('New infection', loc = 'left', y = 0.85, x = 0.02, fontsize = 'medium') |
plt.title('A Daily COVID-19 evolution in MOROCCO') |
x = np.array(df['new_infections'].to_numpy()) |
a = x[100 : 150] |
b = x[151 : 300] |
c = x[301 : 600] |
d = x[600 : 710] |
x_sample = np.array([a[0], b[0], c[0], d[0], x[x.size-1]]) |
print(x_sample) |
y = 1/(1+x_sample*x_sample) |
# print(x) |
poly = lagrange(x_sample, y) |
print(poly) |
Polynomial(poly.coef[::-1]).coef |
# df.plot(x = 'last_updated', y = 'new_infections') |
plt.show() |
x_new = np.arange(-100, 100) / 10 |
plt.scatter(x_sample, y, label = 'data') |
plt.plot(x_new, Polynomial(poly.coef[::-1])(x_new), label = 'Polynomial') |
plt.plot(x_sample, y, label = '1/(1+x*x)') |
plt.legend() |
plt.show() |
Therefore, using Python, we can construct the interpolation polynomial
P(t)=4,04.10−11t4−1,229.10−7t3+8,658.10−5t2−0,01938t+1. | (3.9) |
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 [March3,2020,March3,2020+700day].
Consequently, we inject the function I(t) given by (3.9) (I(t)=P(t)) in Theorem 3 to obtain the function of transmission α(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 [t100,t220], the number of infected people I(t) also increased as the value of α(t) increased, and this period corresponds perfectly to the appearance of the Delta variant in Morroco. Similarly, over the period [t550,t700], 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 α(t) only over a period of time [0,T∗]⊆[0,T] (i.e., T∗<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 α(t) only based on the interpolation polynomial corresponding to a period T∗=300 and n=4 (n is the chosen degree of the interpolator polynomial).
Algorithm 2 An algorithm for obtaining a Lagrange polynomial for any day length and degree |
import pandas as pd |
import json |
import matplotlib.pyplot as plt |
import matplotlib.dates as mdates |
import numpy as np |
from scipy.interpolate import lagrange |
from numpy.polynomial.polynomial import Polynomial |
periods = [300] |
# periods = range(500, 700) |
degrees = [5] |
minimal_degree = 4 |
maximal_degree = 5 |
def split(a, n): |
k, m = divmod(len(a), n) |
return (a[i*k+min(i, m):(i+1)*k+min(i+1, m)] for i in range(n)) |
# Opening JSON file |
f = open('country_newcases.json') |
# returns JSON object as a dictionary |
data = json.load(f) |
df = pd.DataFrame(data) |
# print(df) |
fig, ax = plt.subplots(1, 1, figsize = (10, 7), constrained_layout = True) |
# common to all three: |
ax.plot('last_updated', 'new_infections', data = df) |
# Major ticks every half year, minor ticks every month, |
ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth = (1, 8))) |
ax.xaxis.set_minor_locator(mdates.MonthLocator()) ax.grid(True) |
ax.set_ylabel(r'Cases number') |
for label in ax.get_xticklabels(which = 'major'): |
label.set(rotation = 30, horizontalalignment = 'right') |
ax.ticklabel_format(style = 'plain', useOffset = False, axis = 'y') |
ax.set_title('New infection', loc = 'left', y = 0.85, x = 0.02, fontsize = 'medium') |
plt.title('A Daily COVID-19 evolution in MOROCCO') |
for p in periods: |
for d in degrees: |
x = np.array(df['new_infections'].to_numpy()) |
x_splitted = x[5 : p] |
a = split(x_splitted, d) |
x_sample_list = [] |
for item in list(a): |
x_sample_list.append(item[0]) |
x_sample = np.array(x_sample_list) |
print(x_sample) |
y = 1/(1+x_sample*x_sample) |
poly = lagrange(x_sample, y) |
print(poly) |
In this section, we are interested in the SIR uncontrolled model given by
{dSdt=Λ−αS(t)I(t)−mS(t),dIdt=αS(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t),t∈[0,Tf], | (4.1) |
where S,I and R are the densities of susceptible, infected and removed people respectively, and the transmission coefficient α 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:
{−minimizethenumberofinfectedpeople,−maximizethenumberofremovedpeople, |
and do this with minimal effort and cost. Mathematically, this can be interpreted by the optimization of the function
J(u,v)=∫Tf0[I(t)−R(t)+12ρ1u2(t)+12ρ2v2(t)]dt, | (4.2) |
under the following constraints
{dSdt=Λ−αS(t)I(t)−mS(t)−u(t)S(t),dIdt=αS(t)I(t)−mI(t)−rI(t)−gI(t)−v(t)I(t),dRdt=gI(t)−mR(t)+u(t)S(t)+v(t)I(t).t∈[0,Tf], | (4.3) |
But as mentioned at the beginning of this section the transmission coefficient α 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∗] is devoted to the identifiability of the unknown parameter α (T∗ is as low as possible), and the horizon time [T∗,Tf] 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∈[0,T∗], which has been obtained by interpolation of the database corresponding to T∗ (see Remark 4, Section 3).
Second, we use Theorem 3 to identify the transmission function α(t) corresponding to the following model:
{dSdt=Λ−α(t)S(t)I(t)−mS(t),dIdt=α(t)S(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t).t∈[0,T∗], | (4.4) |
Third, we assign to our unknown parameter α the mean value αm=1T∗∫T∗0α(t)dt, which is calculated by using MATLAB (for the Morrocan case, T∗=60days the Python code gives the prevalence function I(t) which is injected in Theorem 3 to obtain via Matlab the transmission function α(t) for the horizon time [0,T∗] and consequently αm=0.011).
Remark 5. We chose T∗=60, which led us to α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 αm=0.011 is indeed admissible.
T∗ | 40 | 50 | 60 | 70 | 90 | 100 | 110 |
αm | 0.009981 | 0.01504 | 0.0110 | 0.01014 | 0.00993 | 0.0108 | 0.015 |
Having the value of αm=0.011, we use the well identified model,
{dSdt=Λ−αmS(t)I(t)−mS(t)−u(t)S(t),dIdt=αmS(t)I(t)−mI(t)−rI(t)−gI(t)−v(t)I(t),dRdt=gI(t)−mR(t)+u(t)S(t)+v(t)I(t),t∈[T∗,Tf], | (4.5) |
with the initial given data given by
S(T∗)≥0I(T∗)≥0,andR(T∗)≥0, | (4.6) |
to reach our predefined objectives [25], i.e., to minimize the cost function:
J(u,v,Tf)=∫TfT∗[I(t)−R(t)+12ρ1u2(t)+12ρ2v2(t)]dt+ρ3(Tf)2, | (4.7) |
where ρ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 ρ3(Tf)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 Tf as short as possible.
In other words, we search for the optimal controls u∗ and v∗ and the optimal time Tf such that
J(u∗,v∗,Tf)=min{J(u(t),v(t),T)/(u,v)∈UadandT∈R∔}. | (4.8) |
Subject to system (4.5), we take Uad as the control set defined by
Uad={u,v/0<u(t)<umax<∞,0<v(t)<vmax<∞;t∈[T∗,Tf]}. | (4.9) |
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=[f1,f2,f3]⊺ be the right-hand side of the system (4.5) and X=[S,I,R]⊺.
Theorem 6. There exists control functions u∗ and v∗ so that
J(u∗,v∗,Tf)=min(u,v)∈UadJ(u(t),v(t),T). | (4.10) |
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 Uad 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 Uad.
(A5) There exists constants ω1>0,ω2>0 and ψ>1 such that the integrand L(I,R,u,v) of the objective function satisfies
L(I,R,u,v)≥ω2+ω1(∣u∣2+∣v∣2)ψ2. |
(A6) The pay-off term at the optimal time in the objective function φ(X(Tf),Tf)=ρ3(Tf)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 Uad 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≥∣u∣2+∣v∣2 since 0≤u,v≤1.
It is obvious that φ(X(Tf),Tf) is continuous; hence, the last condition is required.
The result follows directly from [26].
It is easy to establish that the corresponding Lagrangian L and Hamiltonian H are, respectively,
L(I,R,u,v,Tf)=I(t)−R(t)+12ρ1u2(t)+12ρ2v2(t)+ρ3T2f, | (4.11) |
H(S,I,R,u,v,λi,Tf)=L(I,R,u,v,Tf)+i=3∑i=1λifi, | (4.12) |
with
{f1=Λ−αmS(t)I(t)−mS(t)−u(t)S(t),f2=αmS(t)I(t)−mI(t)−rI(t)−gI(t)−v(t)I(t),f3=gI(t)−mR(t)+u(t)S(t)+v(t)I(t), |
and where λ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∗,v∗), the optimal time T∗ and the solutions S∗,I∗andR∗ of the corresponding state system (4.5), there exists adjoint variables λ1,λ2 and λ3 satisfying
{˙λ′1=−(−αmI(t)−m−u)λ1−(αmI(t))λ2−uλ3,˙λ′2=−1+αmS(t)λ1−(αmS(t)−m−r−g−v)λ2−(g+v)λ3,λ′3=1+mλ3, | (4.13) |
with the transversality condition λi(Tf)=0 for i=1,2,3.
Furthermore, the optimal control scheme (u∗,v∗) is given by
{u∗(t)=max(min((λ1(t)−λ3(t))ρ1S(t),umax),0),v∗(t)=max(min((λ2(t)−λ3(t))ρ2I(t),vmax),0), | (4.14) |
and the optimal time is given by
T∗f=R(T∗f)−I(T∗f)−12ρ1u(T∗f)2−12ρ2v(T∗f)22ρ3. | (4.15) |
Proof. By applying Pontryagin's maximum principle to the state, we obtain the adjoint equations of the problem by direct computation
{λ′1=−∂H∂S,λ1(Tf)=0,λ′2=−∂H∂I,λ2(Tf)=0,λ′3=−∂H∂R,λ3(Tf)=0, | (4.16) |
and the controls are obtained by calculating the following equations
{∂H∂u=0,∂H∂v=0, | (4.17) |
with the transversality conditions
λi(Tf)=0,i=1,2,3. | (4.18) |
By the bounds in Uad of the controls [30], it is easy to obtain that u∗ and v∗ are given in the form (4.17) of system (4.5).
The transeversality condition for T∗f to be the optimal time can be stated as
H(T∗f,X(T∗f),λ(T∗f),u(T∗f),v(T∗f)=−∂φ(T∗f,X(T∗f)∂t, | (4.19) |
where
φ(T,X(T))=ρ3T2; |
then,
T∗f=R(T∗f)−I(T∗f)−12ρ1u(T∗f)2−12ρ2v(T∗f)22ρ3. | (4.20) |
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∗=60, S(T∗), I(T∗) and R(T∗) 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].
Parameter | Signification | Value |
Λ | Susceptible individual recruitement | 0.525 |
αm | Rate of interaction between susceptibles persons and infected ones | 0.011 |
m | Natural mortality rate | 0.5 |
r | Mortality rate due to the virus | 0.2 |
g | Recovery rate | 0.1 |
S0=S(T∗) | Initial susceptible people | 1000 |
I0=I(T∗) | Initial infected people | 50 |
R0=R(T∗) | Initial recovered people | 15 |
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 Tf=245 days, which was derived from the curve function G(t)=H(t)+∂φ(t)∂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.
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;
{dSdt=Λ−αS(t)I(t)−mS(t),dIdt=αS(t)I(t)−mI(t)−rI(t)−gI(t),dRdt=gI(t)−mR(t).t∈[0,Tf], |
This will be the subject of future work.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Kamal Shah and Thabet Abdeljawad are thankful to Prince Sultan University for paying the APC and for the support through the TAS research lab.
The authors declare that they have no conficts of interest.
[1] |
R. T. Alqahtani, S. Ahmad, A. Akgül, Mathematical analysis of biodegradation model under nonlocal operator in Caputo sense, Mathematics, 9 (2021), 2787. https://doi.org/10.3390/math9212787 doi: 10.3390/math9212787
![]() |
[2] |
M. Arfan, K. Shah, A. Ullah, M. Shutaywi, P. Kumam, Z. Shah, On fractional order model of tumor dynamics with drug interventions under nonlocal fractional derivative, Results Phys., 21 (2021), 103783. https://doi.org/10.1016/j.rinp.2020.103783 doi: 10.1016/j.rinp.2020.103783
![]() |
[3] |
S. Ahmad, A. Ullah, A. Akgül, M. De la Sen, A study of fractional order Ambartsumian equation involving exponential decay kernel, AIMS Math., 6 (2021), 9981–9997. https://doi.org/10.3934/math.2021580 doi: 10.3934/math.2021580
![]() |
[4] | C. Xu, W. Zhang, C. Aouiti, Z. Liu, M. Liao, P. Li, Further investigation on bifurcation and their control of fractional-order bidirectional associative memory neural networks involving four neurons and multiple delays, Math. Method. Appl. Sci., in press, 2022. https://doi.org/10.1002/mma.7581 |
[5] |
C. Xu, Z. Liu, M. Liao, L. Yao, Theoretical analysis and computer simulations of a fractional order bank data model incorporating two unequal time delays, Expert Syst. Appl., 199 (2022), 116859. https://doi.org/10.1016/j.eswa.2022.116859 doi: 10.1016/j.eswa.2022.116859
![]() |
[6] |
C. Xu, W. Zhang, Z. Liu, P. Li, L. Yao, Bifurcation study for fractional-order three-layer neural networks involving four time delays, Cognit. Comput., 14 (2022), 714–732. https://doi.org/10.1007/s12559-021-09939-1 doi: 10.1007/s12559-021-09939-1
![]() |
[7] |
A. Atangana, Fractal-fractional differentiation and integration: Connecting fractal calculus and fractional calculus to predict complex system, Chaos Soliton. Fract., 102 (2017), 396–406. https://doi.org/10.1016/j.chaos.2017.04.027 doi: 10.1016/j.chaos.2017.04.027
![]() |
[8] |
L. Xuan, S. Ahmad, A. Ullah, S. Saifullah, A. Akgül, H. Qu, Bifurcations, stability analysis and complex dynamics of Caputo fractal-fractional cancer model, Chaos Soliton. Fract., 159 (2022), 112113. https://doi.org/10.1016/j.chaos.2022.112113 doi: 10.1016/j.chaos.2022.112113
![]() |
[9] |
S. Ahmad, A. Ullah, A. Akgül, M. De la Sen, Study of HIV disease and its association with immune cells under nonsingular and nonlocal fractal-fractional operator, Complexity, 2021 (2021), 1904067. https://doi.org/10.1155/2021/1904067 doi: 10.1155/2021/1904067
![]() |
[10] |
A. Akgul, N. Ahmed, A. Raza, Z. Iqbal, M. Rafiq, M. A. Rehman, et al., A fractal fractional model for cervical cancer due to human papillomavirus infection, Fractals, 29 (2021), 2140015. https://doi.org/10.1142/S0218348X21400156 doi: 10.1142/S0218348X21400156
![]() |
[11] |
S. Saifullah, A. Ali, K. Shah, C. Promsakon, Investigation of fractal fractional nonlinear Drinfeld–Sokolov–Wilson system with non-singular operators, Results Phys., 33 (2022), 105145. https://doi.org/10.1016/j.rinp.2021.105145 doi: 10.1016/j.rinp.2021.105145
![]() |
[12] |
A. Atangana, M. A. Khan, Fatmawati, Modeling and analysis of competition model of bank data with fractal-fractional Caputo-Fabrizio operator, Alex. Eng. J., 59 (2020), 1985–1998. https://doi.org/10.1016/j.aej.2019.12.032 doi: 10.1016/j.aej.2019.12.032
![]() |
[13] |
S. Ahmad, A. Ullah, A. Akgül, Investigating the complex behaviour of multi-scroll chaotic system with Caputo fractal-fractional operator, Chaos Soliton. Fract., 146 (2021), 110900. https://doi.org/10.1016/j.chaos.2021.110900 doi: 10.1016/j.chaos.2021.110900
![]() |
[14] |
Z. Ahmad, F. Ali, N. Khan, I. Khan, Dynamics of fractal-fractional model of a new chaotic system of integrated circuit with Mittag-Leffler kernel, Chaos Soliton. Fract., 153 (2021), 111602. https://doi.org/10.1016/j.chaos.2021.111602 doi: 10.1016/j.chaos.2021.111602
![]() |
[15] |
S. Qureshi, A. Atangana, A. A. Shaikh, Strange chaotic attractors under fractal-fractional operators using newly proposed numerical methods, Eur. Phys. J. Plus, 134 (2019), 523. https://doi.org/10.1140/epjp/i2019-13003-7 doi: 10.1140/epjp/i2019-13003-7
![]() |
[16] |
K. A. Abro, A. Atangana, Numerical study and chaotic analysis of meminductor and memcapacitor through fractal-fractional differential operator, Arab. J. Sci. Eng., 46 (2021), 857–871. https://doi.org/10.1007/s13369-020-04780-4 doi: 10.1007/s13369-020-04780-4
![]() |
[17] |
Y. Pan, Nonlinear analysis of a four-dimensional fractional hyper-chaotic system based on general Riemann-Liouville-Caputo fractal-fractional derivative, Nonlinear Dyn., 106 (2021), 3615–3636. https://doi.org/10.1007/s11071-021-06951-w doi: 10.1007/s11071-021-06951-w
![]() |
[18] |
S. Saifullah, A. Ali, E. F. D. Goufo, Investigation of complex behaviour of fractal fractional chaotic attractor with mittag-leffler Kernel, Chaos Soliton. Fract., 152 (2021), 111332. https://doi.org/10.1016/j.chaos.2021.111332 doi: 10.1016/j.chaos.2021.111332
![]() |
[19] |
A. Dlamini, E. F. D. Goufo, M. Khumalo, On the Caputo-Fabrizio fractal fractional representation for the Lorenz chaotic system, AIMS Math., 6 (2021), 12395–12421. https://doi.org/10.5465/AMBPP.2021.12421abstract doi: 10.5465/AMBPP.2021.12421abstract
![]() |
[20] |
S. Sampath, S. Vaidyanathan, Ch. K. Volos, V. T. Pham, An eight-term novel four-scroll chaotic system with cubic nonlinearity and its circuit simulation, J. Eng. Sci. Technol., 8 (2015), 1–6. https://doi.org/10.25103/jestr.082.01 doi: 10.25103/jestr.082.01
![]() |
[21] |
J. H. Ma, Y. S. Chen, Study for the bifurcation topological structure and the global complicated character of a kind of nonlinear finance system (I), Appl. Math. Mech., 22 (2001), 1240–1251. https://doi.org/10.1007/BF02437847 doi: 10.1007/BF02437847
![]() |
1. | Zhiliang Li, Lijun Pei, Guangcai Duan, Shuaiyin Chen, A non-autonomous time-delayed SIR model for COVID-19 epidemics prediction in China during the transmission of Omicron variant, 2024, 32, 2688-1594, 2203, 10.3934/era.2024100 | |
2. | Leyla Soudani, Abdelkader Amara, Khaled Zennir, Junaid Ahmad, Duality of fractional derivatives: On a hybrid and non-hybrid inclusion problem, 2024, 0928-0219, 10.1515/jiip-2023-0098 | |
3. | Imane Dehaj, Abdessamad Dehaj, M.A. Aziz-Alaoui, Mostafa Rachik, A new concept of controllability for a class of nonlinear continuous SIR systems, 2025, 192, 09600779, 116013, 10.1016/j.chaos.2025.116013 |
Parameter | Signification |
Λ | Susceptible person recruitement |
α(t) | Disease transmission coefficient |
g | Natural recovery rate |
m | Natural mortality rate |
r | Mortality rate due to the virus |
T∗ | 40 | 50 | 60 | 70 | 90 | 100 | 110 |
αm | 0.009981 | 0.01504 | 0.0110 | 0.01014 | 0.00993 | 0.0108 | 0.015 |
Parameter | Signification | Value |
Λ | Susceptible individual recruitement | 0.525 |
αm | Rate of interaction between susceptibles persons and infected ones | 0.011 |
m | Natural mortality rate | 0.5 |
r | Mortality rate due to the virus | 0.2 |
g | Recovery rate | 0.1 |
S0=S(T∗) | Initial susceptible people | 1000 |
I0=I(T∗) | Initial infected people | 50 |
R0=R(T∗) | Initial recovered people | 15 |
Parameter | Signification |
Λ | Susceptible person recruitement |
α(t) | Disease transmission coefficient |
g | Natural recovery rate |
m | Natural mortality rate |
r | Mortality rate due to the virus |
T∗ | 40 | 50 | 60 | 70 | 90 | 100 | 110 |
αm | 0.009981 | 0.01504 | 0.0110 | 0.01014 | 0.00993 | 0.0108 | 0.015 |
Parameter | Signification | Value |
Λ | Susceptible individual recruitement | 0.525 |
αm | Rate of interaction between susceptibles persons and infected ones | 0.011 |
m | Natural mortality rate | 0.5 |
r | Mortality rate due to the virus | 0.2 |
g | Recovery rate | 0.1 |
S0=S(T∗) | Initial susceptible people | 1000 |
I0=I(T∗) | Initial infected people | 50 |
R0=R(T∗) | Initial recovered people | 15 |