Delay differential equation (DDE) is a type of differential equation in which the derivative of an unknown function at any time depends on the solution at prior times. DDE is mainly used in a large number of system models in physics, biology, and engineering. In addition, DDE is applied in many practical systems such as automatic control, traffic models, neuroscience, and lasers. References [1] and [2] presented a complete study of DDEs and their applications.

On the other hand, when a dynamical phenomenon with a delay is modeled using ordinary DDEs, the model is not always accurate. In general, the initial conditions or parameters of the equations are incomplete or vague. This limitation leads to the study of fuzzy DDEs (FDDEs). In recent years, the theory and application of FDDEs have been developed by many researchers. Lupulescu introduced a DDE with a fuzzy case in [3]. Consequently, Khastan et al. [4] and Hoa et al. [5] presented the results of the existence and uniqueness of solutions of generalized FDDEs. Moreover, Hoa et al. [6] introduced fuzzy delay integro-differential equations with generalized Hukuhara differentiability. Random fuzzy delay integro-differential equations were investigated in [7] and [8].

In many real-world tasks, the memory and hereditary properties of differential operators are broadly perceived to be well anticipated through fractional operators. Fractional order is an extension of the classical order of differentiation and integration. In the last few decades, numerous phenomena in different fields of mathematics, economics, engineering, and science have been more precisely described via fractional derivatives, and fractional differential equations (FDEs) have come up as another powerful device for modeling many difficult types of complex systems. Podlubny [9] provided a major contribution to this field. The existence, uniqueness result, and some integral transforms of FDEs are presented in [9]. Later, Kilbas et al. [10] presented the idea of fractional calculus. A large amount of literature has been dedicated to the theory and applications of FDEs [11–13]. Currently, fractional chaotic systems are gaining considerable interest among many researchers in the different fields of science, mathematics, and engineering. Huang et al. [14] analyzed the behavior of a chaotic system of a fractional order love model with an external environment. Moreover, Huang et al. [15] studied the relationship of Romeo and Juliet love fractional model using a fuzzy function.

The idea of Riemann–Liouville FDE with an uncertainty was introduced by Agarwal et al. [16]. Consequently, many authors presented results on the existence and uniqueness of solutions for fuzzy FDEs (FFDEs) using numerous methods [17–20]. Mazandarani et al. [21] used the modified Euler method to solve FFDEs with a Caputo derivative. Salahshour et al. [22] obtained the solution of FFDEs using the fuzzy Laplace transform method. Ahmed et al. [23] provided analytical and numerical solutions of FFDEs using the Zadeh’s extension principle.

In the present work, our aim is to obtain a solution for FDDEs under a fractional derivative. Many authors solved fuzzy fractional DDEs (FFDDEs) using some numerical methods. For instance, Hoa [24] proposed a solution of FFDDEs that involve the Caputo gH-differentiability using the Adams–Bashforth–Moulton method and also proposed the modified Euler method to solve FFDDEs with an initial condition in [25].

In all the above-cited studies, the source functions and initial conditions are considered as a fuzzy number-valued function. For this function, the derivative is considered as a gH-derivative. This gH-derivative has two types, namely, [(i)-gH] and [(ii)-gH] differentiability. The disadvantages of this derivative are the difficulty in determining which type of derivative to choose and whether the solution is unique. Another difficulty of the gH-derivative is that it may convert the fuzzy problem into two crisp corresponding problems, and we need to choose one solution from a set of solutions that more appropriately fits the problem.

Therefore, our main motivation is to avoid the difficulties in choosing the fuzzy derivative. Then, we enhance the methodology proposed by Fatullayev et al. [26] for FDDEs to FFDDEs. The authors in [26] solved fuzzy problems with an initial condition and a forcing function as a set of fuzzy real functions using the method of Gasilov et al. [27]. Consequently, in [28] the authors provided solutions for fuzzy problems with boundary conditions in which the forcing functions were fuzzy functions.

To the best of our knowledge, this is the first study where this technique is utilized for solving FFDDEs. In the present work, our main objective is to obtain a solution for FDDEs with a fractional derivative. Recently, fuzzy triangular membership function is a beneficial device in many designing applications [29]. Therefore, in this approach, we consider that the initial condition and source functions are in fuzzy triangular forms. Moreover, we prove that if the initial conditions and source functions are fuzzy triangular functions, then the value of the solution is also a fuzzy triangular function. According to this methodology, we present the existence and uniqueness results for of FFDDEs.

This paper is structured as follows: The basic definition and notation of the fractional derivatives, fuzzy sets and fuzzy functions are presented in Section 2. The proposed method for solving FFDDEs is described in Section 3. Numerical applications are performed to test the validity and reliability of the present algorithm in Section 4. Finally, the conclusion is drawn in Section 5.

In this section, we provide certain essential definitions and properties of the fuzzy fractional concept.

Definition 1 [9, 10]

The Riemann–Liouville fractional integral of y(χ) with any order α > 0 is expressed as

Definition 2 [9, 10]

The Caputo derivative of order α (0 < α < 1) of real-valued function y(χ) is defined as

2.1 Fuzzy-Set Theory

We provide fixed set U, and a fuzzy set à on U with membership function μ_à can be defined as à = {(x, μ_Ã(x)) | x ∈ U}.

We also consider a, b, and c as real numbers with a ≤ c ≤ b. The membership function of triangular fuzzy number ũ = (a, c, b) is defined as

According to the geometric explanation of a fuzzy number, c is the vertex, u = a is the left end point, and ū = b is the right end point.

In this study, ũ can be represented by ũ = u_cr + ũ_un, where u_cr is the real part (u_cr = c) and ũ_un is the uncertain part (ũ_un = (a – c, 0, b – c)).

A fuzzy set can also be represented through the definition of r-cut. For each 0 < r ≤ 1, the r-cut of à is expressed as A_r = {x ∈ U | μ_Ã(x) ≥ r}. For r = 0, A₀ = closure(support(Ã)), where support(Ã) = {x ∈ U | μ_Ã(x) > 0}.

The r-cuts of fuzzy number ũ = (a, c, b) are intervals u_r = [u_r, ū_r], where u_r = a + r(c – a) and ū_r = b + r(c – b). This u_r and ū_r can be written as u_r = c + (1 – r)(a – c) and ū_r = c + (1 – r)(b – c). Then, [u_r, ū_r] = c + (1 – r)[a – c, b – c].

2.2 Fuzzy Triangular Function

In our work, we express a fuzzy function as a set of fuzzy real functions in which every real function possesses a specific membership degree. This function is defined by Gasilov et al. [27].

Let F_a(•), F_b(•), and F_c(•) be functions that are continuous on interval I.

Fuzzy set F̃ = 〈F_a, F_c, F_b〉 is called a fuzzy triangular function if it is determined by the following membership function:

The interpretation of triangular fuzzy function F̃ = 〈F_a, F_c, F_b〉 is presented by Gasilov et al in [27].

The fuzzy triangular function is a fuzzy triangular number for each time χ∈I. Then, we can represent it using the following equation:

We consider the following FFDDE:

where m(χ) is a continuous crisp function and τ > 0 is the delay value. Φ̃(χ) = 〈Φ_a, Φ_c, Φ_b〉 and F̃(χ) = 〈F_a, F_c, F_b〉 are triangular fuzzy functions defined on [−τ, 0] and (0, ∞), respectively, and ^CD^α denotes the Caputo derivative operator.

Definition 3 [26]

Fuzzy set Ỹ is the solution of FFDDE (6) on U, and its membership function is defined by

where ϕ(χ) = y(χ).

We can explain Eq. (7) as follows:

Let us consider that y(•) is a function on U. Then, we calculate ϕ(χ) = y(χ) for −τ ≤ χ ≤ 0 and find μ_φ̃(ϕ(•)) = μ₁. Consequently, for given function y(χ), we calculate f(χ) = ^CD^αy(χ) – m(χ)y(χ – τ ) for χ > 0 and find μ_F̃ (f(•)) = μ₂. Next, we determine membership degree μ for function y(•), which is calculated as μ = min{μ₁, μ₂}. Then, we define function set y(•) with membership degree μ as fuzzy solution Ỹ.

As stated in Definition 3, solution Ỹ is a collection of fuzzy real functions y(χ). A function y(χ) has a positive membership degree if it satisfies the following equation:

for ϕ∈supp(Φ̃) and f ∈ support( F̃).

Let us express F̃ as F̃ = f_cr + f̃, where f_cr = F_c, f̃ = 〈f_a, 0, f_b〉 = 〈F_a–F_c, 0, F_b–F_c〉. Let us also express Φ̃ as Φ̃(χ) = ϕ_cr + ϕ, where ϕ_cr = Φ_c, ϕ = 〈ϕ_a, 0, ϕ_b〉 = 〈Φ_a – Φ_c, 0, Φ_b – Φ_c〉.

Because Eq. (6) is linear, we apply the principle of superposition into it.

Then, we split FFDDE (6) into three equations as follows:

• Fractional non-homogeneous equation

  (9){DCαy(χ)=m(χ)y(χ-τ)+fcr(χ),χ>0, 0<α<1,y(χ)=ϕcr(χ),-τ≤χ≤0,

• Fractional homogeneous equation involving an initial function, which is a triangular fuzzy function

  (10){DCαy(χ)=m(χ)y(χ-τ),χ>0, 0<α<1,y(χ)=ϕ˜(χ),-τ≤χ≤0.

• Fractional non-homogeneous equation involving a fuzzy source function and whose initial value is zero

  (11){DCαy(χ)=m(χ)y(χ-τ)+f˜(χ),χ>0, 0<α<1,y(χ)=0,-τ≤χ≤0.

Next, we develop an algorithm to solve these three equations.

3.1 Fractional Non-Homogeneous Equation

We consider fractional non-homogeneous equation, i.e., Eq. (9). We use the method of steps to solve Eq. (9). First, we obtain the solution in the interval [0, τ ]. For 0 ≤ χ ≤ τ, we obtain

By performing fractional integration on both sides of Eq. (12), we obtain unique solution y(χ) = ϕ_cr(0) + J^α[m(χ) ϕ_cr(χ – τ )] + J^α[f_cr(χ)] for 0 ≤ χ ≤ τ, where J^α denotes the fractional integral operators.

3.2 Fractional Homogeneous Equation Involving Initial Fuzzy Function

Next, we solve Eq. (10) in which the initial condition is fuzzy function φ̃ = 〈ϕ_a, 0, ϕ_b〉.

Let us consider that y_a(χ) is the solution of the following equation:

Let us also consider that y_b(χ) is the solution of the following equation: