Positivity and Stability of the Solutions of Caputo Fractional Linear Time-Invariant Systems of Any Order with Internal Point Delays

Abstract

This paper is devoted to the investigation of the nonnegative solutions and the stability and asymptotic properties of the solutions of fractional differential dynamic systems involving delayed dynamics with point delays. The obtained results are independent of the sizes of the delays.

1. Introduction

The theory of fractional calculus is basically concerned with the calculus of integrals and derivatives of any arbitrary real or complex orders. In this sense, it may be considered as a generalization of classical calculus which is included in the theory as a particular case. The former ideas have been stated about three hundred years ago, but the main mathematical developments and applications of fractional calculus have been of increasing interest from the seventies. There is a good compendium of the state of the art of the subject and the main related existing mathematical results with examples and case studies in [1]. There are a lot of results concerning the exact and approximate solutions of fractional differential equations of Riemann-Liouville and Caputo types, [1–4], fractional derivatives involving products of polynomials, [5, 6], fractional derivatives and fractional powers of operators, [7–9], boundary value problems concerning fractional calculus (see, e.g., [1, 10]), and so forth. There is also an increasing interest in the recent mathematical literature in the characterization of dynamic fractional differential systems oriented towards several fields of science like physics, chemistry or control theory because it is a powerful tool for later applications in all fields requiring support via ordinary, partial derivatives, and functional differential equations. Perhaps the reason of interest of fractional calculus is that the numerical value of the fraction parameter allows a closer characterization of eventual uncertainties present in the dynamic model compared to the alternative use of structured uncertainties. We can find, in particular, a lot of literature concerned with the development of Lagrangian and Hamiltonian formulations where the motion integrals are calculated though fractional calculus and also in related investigations concerning dynamic and damped and diffusive systems [11–17] as well as the characterization of impulsive responses or its use in applied optics related, for instance, to the formalism of fractional derivative Fourier plane filters (see, e.g., [16–18]) and Finance [19]. Fractional calculus is also of interest in control theory concerning, for instance, heat transfer, lossless transmission lines, the use of discretizing devices supported by fractional calculus, and so forth (see, e.g., [20–22]). In particular, there are several recent applications of fractional calculus in the fields of filter design, circuit theory and robotics, [21, 22], and signal processing, [17]. Fortunately, there is an increasing mathematical literature, currently available on fractional differ-integral calculus, which can formally support successfully the investigations in other related disciplines.

This paper is concerned with the investigation of the solutions of time-invariant fractional differential dynamic systems, [23, 24], involving point delays which leads to a formalism of a class of functional differential equations, [25–31]. Functional equations involving point delays are a crucial mathematical tool to investigate real process where delays appear in a natural way like, for instance, transportation problems, war and peace problems, or biological and medical processes. The main interest of this paper is concerned with the positivity and stability of solutions independent of the sizes of the internal delays and also with obtaining results being independent of the eventual mutual coincidence of some values of delays, [31–33]. It has to be pointed out that the positivity of the solutions is a crucial property in investigating some dynamic systems like biological systems or epidemic models, [32, 33], where positivity is an essential requirement since negative solutions have nonsense at any time instant. It is also a relevant property concerning the existence and characterization of oscillatory solutions of differential equations, [34]. Most of the results are centred in characterizations via Caputo fractional differentiation although some extensions are presented concerned with the classical Riemann-Liouville differ integration. It is proved that the existence of nonnegative solutions independent of the sizes of the delays and the stability properties of linear time-invariant fractional dynamic differential systems subject to point delays may be characterized with sets of precise mathematical results.

1.1. Notation

𝐙, 𝐑, and 𝐂 are the sets of integer, real, and complex numbers, 𝐙₊ and 𝐑₊ are the positive integer and real numbers, and 𝐙₀₊∶=𝐙₊𝐂∪{0},𝐑₀₊∶=𝐑₊∪{0},𝐂₊∶={𝑧∈𝐂∶Re𝑧>0},₀₊∶={𝑧∈𝐂∶Re𝑧≥0},𝑛∶={1,2,…,𝑛}.(1.1)

The following notation is used to characterize different levels of positivity of matrices: 𝐑^𝑛×𝑚0+∶={𝑀=(𝑀_𝑖𝑗)∈𝐑^𝑛×𝑚∶𝑀_𝑖𝑗≥0;∀(𝑖,𝑗)∈𝑛×𝑚} is the set of all 𝑛×𝑚 real matrices of nonnegative entries. If 𝑀∈𝐑^𝑛×𝑚 then 𝑀≥0 is used as a simpler notation for 𝑀∈𝐑^𝑛×𝑚0+.

𝐑^+𝑛×𝑚∶={0≠𝑀=(𝑀_𝑖𝑗)∈𝐑^𝑛×𝑚∶𝑀_𝑖𝑗≥0;∀(𝑖,𝑗)∈𝑛×𝑚} is the set of all nonzero 𝑛×𝑚 real matrices of nonnegative entries (i.e., at least one of their entries is positive). If 𝑀∈𝐑^𝑛×𝑚 then 𝑀>0 is used as a simpler notation for 𝑀∈𝐑^+𝑛×𝑚.

𝐑^𝑛×𝑚++∶={𝑀=(𝑀_𝑖𝑗)∈𝐑^𝑛×𝑚∶𝑀_𝑖𝑗>0;∀(𝑖,𝑗)∈𝑛×𝑚} is the set of all 𝑛×𝑚 real matrices of positive entries. If 𝑀∈𝐑^𝑛×𝑚 then 𝑀≫0 is used as a simpler notation for 𝑀∈𝐑^𝑛×𝑚++. The superscript 𝑇 denotes the transpose, 𝑀^𝑇𝑖 and 𝑀_𝑗 are, respectively, the 𝑖th row and the 𝑗th column of the matrix 𝑀.

A close notation to characterize the positivity of vectors is the following: 𝐑^𝑛0+∶={𝑣=(𝑣₁,𝑣₂,…,𝑣_𝑛)^𝑇∈𝐑^𝑛∶𝑣_𝑖≥0;∀𝑖∈𝑛} is the set of all 𝑛 real vectors of nonnegative components. If 𝑣∈𝐑^𝑛 then 𝑣≥0 is used as a simpler notation for 𝑣∈𝐑^𝑛0+. 𝐑^𝑛+∶={0≠𝑣=(𝑣₁,𝑣₂,…,𝑣_𝑛)^𝑇∈𝐑^𝑛∶𝑣_𝑖≥0;∀𝑖∈𝑛} is the set of all 𝑛 real nonzero vectors of nonnegative components (i.e., at least one component is positive). If 𝑣∈𝐑^𝑛 then 𝑣>0 is used as a simpler notation for 𝑣∈𝐑^𝑛+.

𝐑^𝑛++∶={𝑣=(𝑣₁,𝑣₂,…,𝑣_𝑛)^𝑇∈𝐑^𝑛∶𝑣_𝑖>0;∀𝑖∈𝑛} is the set of all 𝑛 real vectors of positive components. If 𝑣∈𝐑^𝑛 then 𝑣≫0 is used as a simpler notation for 𝑣∈𝐑^𝑛++.

𝑀=(𝑀_𝑖𝑗)∈𝐑^𝑛×𝑛 is a Metzler matrix if 𝑀_𝑖𝑗≥0; for all (𝑖,𝑗≠𝑖)∈𝑛×𝑛. 𝑀𝐑^𝑛×𝑛 is the set of Metzler matrices of order 𝑛.

The maximum real eigenvalue, if any, of a real matrix 𝑀, is denoted by 𝜆_max(𝑀). Multiple subscripts of vector, matrices, and vector and matrix functions are separated by commas only in the case that, otherwise, some confusion could arise as, for instance, when some of the subscripts is an expression involving several indices.

2. Some Background on Fractional Differential Systems

Assume that 𝑓∶[𝑎,𝑏]→𝐂^𝑛 for some real interval [𝑎,𝑏]⊂𝐑 satisfies𝑓∈𝐶^𝑘−2((𝑎,𝑏),𝐑^𝑛) and, furthermore, 𝑑^𝑘−1𝑓(𝜏)/𝑑𝜏^𝑘−1 exists everywhere in [𝑎,𝑏] for 𝑘=[Re𝛼]+1 for some 𝛼∈𝐂₀₊. Then, the Riemann-Liouville left-sided fractional derivative ^RL𝐷^𝛼𝑎+𝑓 of order 𝛼∈𝐂₀₊ of the vector function 𝑓 in [𝑎,𝑏] is pointwise defined in terms of the Riemann-Liouville integral as^RL𝐷^𝛼𝑎+𝑓(𝑡)∶=Γ(𝑘1−𝛼)𝑑^𝑘𝑑𝑡^𝑘^𝑡𝑎𝑓(𝜏)(𝑡−𝜏)^{𝛼+1−𝑘}𝑑𝜏,𝑡∈[𝑎,𝑏],(2.1) where the integer 𝑘 is given by 𝑘=[Re𝛼]+1 and Γ∶𝐂⧵𝐙₀₋→𝐂, where 𝐙₀₋∶={𝑛∈𝐙∶𝑛≤0}, is the Γ-function defined by Γ(𝑧)∶=∫^∞0𝜏^𝑧−1𝑒^−𝜏𝑑𝜏; 𝑧∈𝐂⧵𝐙₀₋. If 𝑓∈𝐶^𝑘−1((𝑎,𝑏),𝐑^𝑛) and, furthermore, 𝑓^(𝑘)(𝜏)≡𝑑^𝑘𝑓(𝜏)/𝑑𝜏^𝑘 exists everywhere in [𝑎,𝑏], then the Caputo left-sided fractional derivative ^𝐶𝐷^𝛼𝑎+𝑓 of order 𝛼∈𝐂₀₊ of the vector function 𝑓 in [𝑎,𝑏] is pointwise defined in terms of the Riemann-Liouville integral as^𝐶𝐷^𝛼𝑎+𝑓(𝑡)∶=Γ(𝑘1−𝛼)^𝑡𝑎𝑓^(𝑘)(𝜏)(𝑡−𝜏)^{𝛼+1−𝑘}𝑑𝜏,𝑡∈[𝑎,𝑏],(2.2) where 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₀₊ and 𝑘=𝛼 if 𝛼∈𝐙₀₊. The following relationship between both fractional derivatives holds provided that they exist (i.e., if 𝑓∶[𝑎,𝑏]→𝐂^𝑛 possesses Caputo left-sided fractional derivative in [𝑎,𝑏]), [1]^𝐶𝐷^𝛼𝑎+𝑓(𝑡)=^RL𝐷^𝛼𝑎+𝑓(𝜏)−^𝑘−1_𝑗=0𝑓^(𝑗)(𝑎)(𝜏−𝑎)^𝑗𝑗!(𝑡),𝑡∈[𝑎,𝑏].(2.3) Since Re𝛼≤𝑘, the above formula relating both fractional derivatives proves the existence of the Caputo left-sided fractional derivative in [𝑎,𝑏] if the Riemann-Liouville one exists in [𝑎,𝑏].

3. Solution of a Fractional Differential Dynamic System of Any Order 𝛼 with Internal Point Delays

Consider the linear and time-invariant differential functional Caputo fractional differential system of order 𝛼:^𝐶𝐷^𝛼0+𝑥(𝑡)=Γ(𝑘1−𝛼)^𝑡𝑎𝑓^(𝑘)(𝜏)(𝑡−𝜏)^{𝛼+1−𝑘}𝑑𝜏=^𝑝_𝑖=0𝐴_𝑖𝑥𝑡−ℎ_𝑖+𝐵𝑢(𝑡)(3.1) with 𝑘−1<𝛼(∈𝐑₊)≤𝑘; 𝑘−1, 𝑘∈𝐙₀₊, 0=ℎ₀<ℎ₁<ℎ₂<⋯<ℎ_𝑝=ℎ<∞ being distinct constant delays, 𝐴₀,𝐴_𝑖∈𝐑^𝑛×𝑛(𝑖∈𝑝∶={1,2,…,𝑝}), are the matrices of dynamics for each delay ℎ_𝑖,𝑖∈𝑝∪{0}, 𝐵∈𝐑^𝑛×𝑚 is the control matrix. The initial condition is given by 𝑘  𝑛-real vector functions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛, with 𝑗∈𝑘−1∪{0}, which are absolutely continuous except eventually in a set of zero measure of [−ℎ,0]⊂𝐑 of bounded discontinuities with 𝜑_𝑗(0)=𝑥_𝑗(0)=𝑥^(𝑗)(0)=𝑥_𝑗0. The function vector 𝑢∶𝐑₀₊→𝐑^𝑚 is any given bounded piecewise continuous control function. The following result is concerned with the unique solution on 𝐑₀₊ of the above differential fractional system (3.1). The proof follows directly from a parallel existing result from the background literature on fractional differential systems by grouping all the additive forcing terms of (3.1) in a unique one (see, e.g., [1, (1.8.17), (3.1.34)–(3.1.49)], with 𝑓(𝑡)≡∑^𝑝𝑖=1𝐴_𝑖𝑥(𝑡−ℎ_𝑖)+𝐵𝑢(𝑡)).

Theorem 3.1. The linear and time-invariant differential functional fractional differential system (3.1) of any order 𝛼∈𝐂₀₊ has a unique solution on 𝐑₀₊ for each given set of initial functions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛, 𝑗∈𝑘−1∪{0} being absolutely continuous except eventually in a set of zero measures of [−ℎ,0]⊂𝐑 of bounded discontinuities with 𝜑_𝑗(0)=𝑥_𝑗(0)=𝑥^(𝑗)(0)=𝑥_𝑗0; 𝑗∈𝑘−1∪{0} and each given control 𝑢∶𝐑₀₊→𝐑^𝑚 being a bounded piecewise continuous control function. Such a solution is given by 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Φ^𝛼𝑗0(𝑡)𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀+Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏_𝑝_𝑖=1^𝑡ℎ_𝑖Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝑥_𝛼𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Φ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏,𝑡∈𝐑₀₊(3.2) with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼 if 𝛼∈𝐙₊, Φ_𝛼𝑗0𝐸(𝑡)∶=𝑡^𝑗𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼,Φ_𝛼(𝑡)∶=𝑡^𝛼−1𝐸_𝛼,𝛼𝐴₀𝑡^𝛼,(3.3)_𝛼,𝑗𝐴₀𝑡^𝛼∶=^∞_ℓ=0𝐴₀𝑡^𝛼^ℓΓ(𝛼ℓ+𝑗),𝑗∈𝑘−1∪{0,𝛼}(3.4) for 𝑡≥0 and Φ_𝛼0(𝑡)=Φ_𝛼(𝑡)=0 for 𝑡<0, where 𝐸_𝛼,𝑗(𝐴₀𝑡^𝛼) are the Mittag-Leffler functions.

Now consider that the right-hand side of (3.1) is the evaluation of a Riemann-Liouville fractional differential system of the same order 𝛼 as follows: ^RL𝐷^𝛼0+𝑥(𝑡)=^𝑝_𝑖=0𝐴_𝑖𝑥𝑡−ℎ_𝑖+𝐵𝑢(𝑡)(3.5) under the same functions of initial conditions as those of (3.1). Through the formula (2.3) relating Caputo and Riemann-Liouville left-sided fractional derivatives of the same order 𝛼, one gets^𝐶𝐷^𝛼0+𝑥(𝑡)=^𝑝_𝑖=0𝐴_𝑖𝑥𝑡−ℎ_𝑖+𝐵𝑢(𝑡)−^RL𝐷^𝛼0+^𝑘−1_𝑗=0𝑥_𝑗0𝜏^𝑗𝑗!(𝑡).(3.6) Since the Caputo left-sided fractional derivative and the Riemann-Liouville fractional integral of order 𝛼∈𝐂₊∶={𝐙₊∪{𝑧∈𝐂₊∶Re𝑧∉𝐙₊}} are inverse operators (what is not the case if 𝛼∉𝐂₊), (see [1, Lemma 2.21(a)]), one gets from (3.6), (2.3), and (3.2) if 𝛼∈𝐂 the subsequent result for the fractional differential system (3.5) on 𝐑₀₊.

Corollary 3.2. If (3.5) of any order 𝛼∈𝐂₊ is replaced with (3.1) under the same initial conditions then its unique solution on 𝐑₀₊ is given by 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Φ^𝛼𝑗0(𝑡)−𝑡^𝑗𝑗!𝐼^𝑛𝑥^𝑗0+^𝑝_𝑖=1^ℎ𝑖₀+Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏_𝑝_𝑖=1^𝑡ℎ_𝑖Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝑥_𝛼𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Φ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏,𝛼∉𝐙₀₊;𝑡∈𝐑₀₊=^𝑘−1_𝑗=0𝐸^𝛼,𝑗+1(𝑡)−1𝑗!𝐼^𝑛𝑡^𝑗𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀+Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏_𝑝_𝑖=1^𝑡ℎ_𝑖Φ_𝛼(𝑡−𝜏)𝐴_𝑖𝑥_𝛼𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Φ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏,𝛼∉𝐙₀₊,𝑡∈𝐑₀₊(3.7) with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼 if 𝛼∈𝐙₊.

Another mild evolution operator can be considered to construct the unique solution of (3.1) by considering the control effort as the unique forcing term of (3.1) and the functions of initial conditions as forcing terms. See the corresponding expressions obtainable from [1, (1.8.17), (3.1.34)–(3.1.49)], with the identity (𝑓(𝑡)≡𝐵𝑢(𝑡)) and the evolution operator defined in [2, 3] for the standard (nonfractional differential system), that is, 𝛼=1 in (3.1). Thus, another equivalent expression for the unique solution of the Caputo fractional differential system of order 𝛼 is given in the subsequent result.

Theorem 3.3. The solution of (3.1) given in Theorem 3.1 is equivalently rewritten as follows: 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Ψ^𝛼𝑗0(𝑡)𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀Ψ_𝛼𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Ψ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏(3.8) for 𝑡∈𝐑₀₊, any 𝛼∈𝐂₊ with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼 if 𝛼∈𝐙₊; Ψ_𝛼𝑗0(𝑡)∶=𝑡^𝑗𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼+^𝑝_𝑖=1^𝑡0Ψ𝜏^𝛼−1𝐸_𝛼𝛼𝐴₀𝜏^𝛼𝐴_𝑖Ψ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖𝑑𝜏,_𝛼(𝑡)∶=𝑡^𝛼−1𝐸_𝛼𝛼𝐴₀𝑡^𝛼+^𝑝_𝑖=1^𝑡0𝜏^𝛼−1𝐸_𝛼𝛼𝐴₀𝜏^𝛼𝐴_𝑖Ψ_𝛼𝑡−𝜏−ℎ_𝑖𝑑𝜏(3.9) for 𝑡≥0 and Ψ_𝛼𝑗0(𝑡)=Ψ_𝛼(𝑡)=0, 𝑗∈𝑘−1∪{0} for 𝑡∈[−ℎ,0).

Also, the solution to the Riemann-Liouville fractional differential system (3.5) under the same initial conditions as those of (3.4) is given in the next result for 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ based on (3.6).

Corollary 3.4. If (3.5) being of order 𝛼∈𝐂₊ is replaced with (3.1) under the same initial conditions then its unique solution on 𝐑₀₊ is given by 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Ψ^𝛼𝑗0(𝑡)−𝑡^𝑗𝑗!𝐼^𝑛𝑥^𝑗0+^𝑝_𝑖=1^ℎ𝑖₀Ψ_𝛼𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Ψ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏(3.10) with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼 if 𝛼∈𝐙₊ which is identical to that given in Corollary 3.2.

Particular cases of interest of the solution of (3.1) given in Theorem 3.3 are(1)𝛼=𝑘 which yields the solution:𝑥_𝑘(𝑡)=^𝑘−1_𝑗=0Ψ^𝑘𝑗0(𝑡)𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀Ψ_𝑘𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Ψ_𝑘(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏,(3.11)(2)a further particular case 𝛼=𝑘=1 which yields the solution: 𝑥₁(𝑡)=Ψ₁(𝑡)𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀Ψ₁(𝑡−𝜏)𝜑₀𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Ψ₁(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏(3.12) since Ψ₁₀₀(𝑡)=Ψ₁(𝑡), 𝑡∈𝐑₀₊ which is the unique solution of (𝐷𝑥)(𝑡)=∑^𝑝𝑖=0𝐴_𝑖(𝑡−ℎ_𝑖)+𝐵𝑢(𝑡) under any almost everywhere absolutely continuous function (except eventually in some subset of zero measure of [−ℎ,0] of bounded discontinuities) of initial conditions  𝜑≡𝜑₀∶[−ℎ,0]→𝐑^𝑛. Use for this case, the less involved notations (Ψ𝑥)(𝑡)=(Ψ₁₀₀𝑥)(𝑡)=(Ψ₁𝑥)(𝑡) for the smooth evolution operator from 𝐑₀₊×𝐑^𝑛 to 𝐑^𝑛, and Φ(𝑡)=Φ₁₀₀(𝑡)=Φ₁(𝑡)=𝑒^𝐴0𝑡, 𝑡∈𝐑₊ for the exponential matrix function 𝑒^𝐴0𝑡 from 𝐑₀₊ to 𝐑^𝑛×𝑛, which defines a 𝐶₀-semigroup (𝑒^𝐴0𝑡,𝑡∈𝐑₀₊) of infinitesimal generator 𝐴₀ from 𝐑₀₊ to 𝐿(𝐑^𝑛). Then, the unique solution 𝑥(𝑡)≡𝑥₁(𝑡),𝑡∈𝐑₊ for the given function of initial conditions is𝑥(𝑡)=Φ(𝑡)𝑥₀+^𝑝_𝑖=1^ℎ𝑖₀+_𝑝Φ(𝑡−𝜏)𝐴_𝑖𝜑𝜏−ℎ_𝑖_𝑖=1^𝑡ℎ_𝑖=Ψ(𝑡)𝑥₀Φ(𝑡−𝜏)𝐴_𝑖𝑥𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Φ(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏+^𝑝_𝑖=1^ℎ𝑖₀Ψ(𝑡−𝜏)𝜑𝜏−ℎ_𝑖𝑑𝜏+^𝑡0Ψ(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏,𝑡∈𝐑₀₊(3.13) and 𝑥(𝑡)=𝜑(𝑡) for 𝑡∈[−ℎ,0], where Φ(𝑡)=𝑒^𝐴0𝑡 satisfies Φ̇(𝑡)=𝐴₀Φ(𝑡)𝑡∈𝐑, and Ψ̇(𝑡)=∑^𝑝𝑖=0𝐴_𝑖Ψ(𝑡−ℎ_𝑖) with Ψ(0)=Φ(0)=𝐼_𝑛 (the 𝑛-identity matrix) and Ψ(𝑡)=0, 𝑡∈[−ℎ,0) which has a unique solution Ψ(𝑡)=𝑒^𝐴0𝑡(𝐼_𝑛+∑^𝑝𝑖=1∫^𝑡ℎ_𝑖𝑒^−𝐴0𝜏𝐴_𝑖Ψ(𝜏−ℎ_𝑖)𝑑𝜏) for 𝑡∈𝐑₀₊, [2, 3, 25, 26].

A problem of interest when considering a set of 𝑝 delays in [0,ℎ] is the case of potentially repeated delays, then subject to 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝≤ℎ<∞, with 𝑞 of them {ℎ_𝑗𝑝,𝑗∈𝑞∪{0}} being distinct, each being repeated 1≤𝜈_𝑗≤𝑝(𝑗∈𝑞∪{0}) times so that0=ℎ₀=ℎ_0𝑝<ℎ_1𝑝<⋯<ℎ_𝑞𝑝≤ℎ<∞,^𝑝_𝑗=0ℎ𝜈_𝑗𝑝=𝑝+1,_𝑗𝑝=ℎ_𝑘+𝑖,𝑘=^𝑗−1_ℓ=0𝜈_ℓ,∀𝑖∈𝜈_𝑗,𝑗∈𝑞∪{0}.(3.14)

Thus, the following result holds from Theorem 3.3 by grouping the terms of the delayed dynamics corresponding to the same potentially repeated delays.

Theorem 3.5. The Caputo solutions to the subsequent Caputo and Riemann-Liouville fractional differential systems of order 𝛼 with 𝑝≥0 (potentially repeated) delays and 0≤𝑞≤𝑝 distinct delays: ^𝐶𝐷^𝛼0+𝑥(𝑡)=^𝑝_𝑖=0𝐴_𝑖𝑡−ℎ_𝑖+𝐵𝑢(𝑡),^RL𝐷^𝛼0+𝑥(𝑡)=^𝑝_𝑖=0𝐴_𝑖𝑡−ℎ_𝑖+𝐵𝑢(𝑡)(3.15) on 𝐑₀₊ for the given set of initial conditions on [−ℎ,0) are given by 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Ψ^𝛼𝑗0(𝑡)𝑥_𝑗0+^𝑞_𝑖=1^ℎ𝑖𝑝₀Ψ_𝛼𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑝𝑑𝜏+^𝑡0Ψ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏(3.16) for any 𝛼∈𝐂₊ with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₀₊ and 𝑘=𝛼 if 𝛼∈𝐙₀₊, and, respectively by 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Ψ^𝛼𝑗0(𝑡)−𝑡^𝑗𝑗!𝑥^𝑗0+^𝑞_𝑖=1^ℎ𝑖𝑝₀Ψ_𝛼𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑝𝑑𝜏+^𝑡0Ψ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏(3.17) for any 𝛼∈𝐂₊ with 𝑘=[Re𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼 if 𝛼∈𝐙₊, where Ψ_𝛼𝑗0(𝑡)∶=𝐸_𝛼,𝑗+1⎛⎜⎜⎝⎛⎜⎜⎝^𝜈0−1_𝑖=0𝐴_𝑖+_𝑞⎞⎟⎟⎠𝑡_𝛼⎞⎟⎟⎠_𝑖=1^𝑡0𝜏^𝛼−1𝐸_𝛼,𝛼⎛⎜⎜⎝⎛⎜⎜⎝^𝜈0−1_𝑖=0𝐴_𝑖⎞⎟⎟⎠𝜏_𝛼⎞⎟⎟⎠^𝜈𝑖_ℓ=1𝐴_(∑^{𝑖−1𝑗=0}_{𝜈𝑗+ℓ)}Ψ^𝛼𝑗0𝑡−𝜏−ℎ_𝑖𝑑𝜏,(3.18)𝑡≥0 and Ψ_𝛼𝑗0(𝑡)=Ψ_𝛼(𝑡)=0, 𝑗∈𝑘−1∪{0} for 𝑡∈[−ℎ,0).

4. Nonnegativity of the Solutions

The positivity of the solutions of (3.1) independent of the values of the delays is now investigated under initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0}.

Theorem 4.1. The Caputo fractional differential system (3.1) under the delay constraint 0=ℎ₀<ℎ₁<ℎ₂<⋯<ℎ_𝑝=ℎ<∞ for any given absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0} and any piecewise continuous vector function 𝑢∶𝐑₀₊→𝐑^𝑛0+ if 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊,forall𝑡∈𝐑₀₊; for all 𝛼∈𝐑₊ has following properties:
(i)Φ_𝛼𝑗0(𝑡) is nonsingular; forall𝑗∈𝑘−1∪{0} and Φ_𝛼(𝑡)≥0; for all 𝑡∈𝐑₀₊ (if 𝐵∈𝐑^+𝑛×𝑚 then Φ_𝛼(𝑡)>0; for all 𝑡∈𝐑₀₊),(ii) (1)𝐴₀∈𝑀𝐑^𝑛×𝑛⇔Φ(𝑡)≡Φ₁₀₀(𝑡)>0; for all 𝑡∈𝐑₀₊,(2)𝐴₀∈𝑀𝐑^𝑛×𝑛⇒Φ_𝛼𝑗0(𝑡)>0;forall𝑗∈𝑘−1∪{0}; for all 𝑡∈[0,𝑡_𝑗) for some sufficiently small 𝑡_𝑗∈𝐑₊ with Φ_𝛼00(𝑡)>0, for all 𝑡∈𝐑₀₊(i.e.,𝑡₀=∞). This property holds for all 𝑡∈𝐑₀₊  (i.e., 𝑡_𝑗=∞; for all 𝑗∈𝑘−1∪{0}) if, in addition, either 𝐴₀≥0 or if 𝐴₀ is nilpotent or if 0<𝛼≤𝑘=1. Furthermore, there are at least n entries (one per row) of Φ_𝛼𝑗0(𝑡) being positive; forall𝑡∈𝐑₀₊; (iii)Any solution (3.2) to any Caputo fractional differential system (3.1) is nonnegative independent of the delays; that is, 𝑥_𝛼(𝑡)∈𝐑^𝑛0+; for all 𝑡∈[−ℎ,𝑡)∩𝐑₀₊ for some 𝑡∈𝐑₀₊, for any set of delays satisfying 0=ℎ₀<ℎ₁<ℎ₂<⋯<ℎ_𝑝≤ℎ<∞ and any absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, for all 𝑗∈𝑘−1∪{0} and any piecewise continuous control 𝑢∶𝐑₀₊→𝐑^𝑚0+, if and only if 𝐴₀∈𝑀𝐑^𝑛×𝑛 for 𝑡∈𝐑₀₊ being sufficiently small. Furthermore, 𝑥_𝛼(𝑡)∈𝐑^𝑛0+; for all 𝑡∈[−ℎ,0)∪𝐑₀₊ if, in addition, either 𝐴₀≥0 or if 𝐴₀ is nilpotent or if 0<𝛼≤𝑘=1,𝐴_𝑖∈𝐑^𝑛×𝑛0+(∀𝑖∈𝑝) and 𝐵∈𝐑^𝑛×𝑚0+.

Proof. It is now proven that Φ_𝛼𝑗0(𝑡)≥0; for all 𝑡∈𝐑₀₊⇒𝐴₀∈𝑀𝐑^𝑛×𝑛; for all 𝛼∈𝐑₊ for any 𝑗∈𝑘−1∪{0}. First, note the following. If 𝛼=𝑘=1 then Φ_𝛼00(𝑡)=𝐸_1,1(𝐴₀𝑡)=Φ(𝑡)=∑^∞ℓ=0𝐴^ℓ0𝑡^ℓ/ℓ!=𝑒^𝐴0𝑡≥0 if 𝐴₀∈𝑀𝐑^𝑛×𝑛 from the above part of the proof and also 𝐴₀∈𝑀𝐑^𝑛×𝑛⇒Φ(𝑡)≥0; for all 𝑡∈𝐑₀₊. This follows by contradiction. Assume that Φ_𝑖𝑚(𝑡)<0 for some 𝑡∈𝐑₊. Consider the positive differential system 𝑥̇(𝑡)=𝐴₀𝑥(𝑡), 𝑥(0)=𝑒_𝑗, 𝐴₀∈𝑀𝐑^𝑛×𝑛 so that 𝑥_𝑖(𝑡)=−|Φ_𝑖𝑚(𝑡)|<0 which contradicts the system being positive. Thus, 𝐴₀∈𝑀𝐑^𝑛×𝑛⇔Φ(𝑡)≥0; for all 𝑡∈𝐑₀₊. Furthermore, since Φ(𝑡) is a fundamental matrix of solutions of the differential system, it is non-singular for all finite time and the above result is weakened as follows:
𝐴₀∈𝑀𝐑^𝑛×𝑛⇔(Φ(𝑡)=Φ_𝛼00(𝑡)=𝑒^𝐴0𝑡>0∧Φ(𝑡) is non-singular; for all 𝑡∈𝐑₀₊). Since Φ(𝑡) is nonsingular; for all 𝑡∈𝐑₀₊ at least 𝑛 of its entries (one per-row) is positive. Property (i) has been proven. Now, one gets from (3.3)-(3.4): Φ_𝛼𝑗0(𝑡)=𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼=^∞_ℓ=0𝐴^ℓ0𝑡^𝛼ℓΓ(𝛼ℓ+𝑗+1)=^∞_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑡^{(𝛼−1)ℓ}ℓ!(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑗),𝑗∈𝑘−1∪{0}.(4.1) Let 𝑒_𝑖 the 𝑖th unit Euclidean vector of 𝐑^𝑛 whose 𝑖th component is 1. Then, one obtains for all 𝑗∈𝑘−1∪{0}, irrespective of the value of 𝛼∈𝐑₊ and 𝑘∈𝐙₊ being 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, provided that 𝐴₀∈𝑀𝐑^𝑛×𝑛: Φ_{𝛼𝑗0𝑖𝑚}(𝑡)=𝑒^𝑇𝑖Φ_𝛼𝑗0(𝑡)𝑒_𝑚=𝑒^𝑇𝑖^∞_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑡^{(𝛼−1)ℓ}ℓ!Γ(𝛼ℓ+𝑗+1)𝑒^𝑚=𝑒^𝑇𝑖^∞_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑡^{(𝛼−1)ℓ}ℓ!(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑗)𝑒^𝑚≥𝑒^𝑇𝑖𝑒^𝐴0𝑡𝑒_𝑚min_{0≤ℓ≤∞}𝑡^{(𝛼−1)ℓ}ℓ!=𝑒_𝑇𝑖Γ(𝛼ℓ+𝑘)^∞_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑡^{(𝛼−1)ℓ}ℓ!(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑗)𝑒^𝑚(4.2)=𝑒^𝑇𝑖^𝑁_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑒^𝑚min_{0≤ℓ≤𝑁}𝑡^{(𝛼−1)ℓ}ℓ!=𝑒_𝑇𝑖(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑘)𝑒^𝐴0𝑡𝑒_𝑚min_{0≤ℓ≤𝑁}𝑡^{(𝛼−1)ℓ}Γ(ℓ+1)≥𝑒_𝑇𝑖(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑘)𝑒^𝐴0𝑡𝑒_𝑚min_{0≤ℓ≤𝑁}𝑡^{(𝛼−1)ℓ}Γ(ℓ+1)(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑘)(4.3) for all (𝑖,𝑚)∈𝑛×𝑛, for all 𝑡∈𝐑₀₊ since 𝐴₀∈𝑀𝐑^𝑛×𝑛⇔Φ(𝑡)≥0; for all 𝑡∈𝐑₀₊, for some 𝑁(≤∞)∈𝐙₀₊ and 𝑁 is finite if and only if 𝐴₀ is nilpotent (of degree 𝑁). Equation (4.3) implies that 𝑒^𝑇𝑖(𝑒^𝐴0𝑡)𝑒_𝑚≥0 and then (Φ_{𝛼𝑗0𝑖𝑚}(𝑡))≥0, for all (𝑖,𝑚)∈𝑛×𝑛 in the following cases:
(a) 𝑁≤∞, 𝑡∈[0,𝑡), since ∩𝐑₀₊ and some sufficiently small 𝑡∈𝐑₊, since 𝐴₀∈𝑀𝐑^𝑛×𝑛⇒𝐼+𝐴₀𝑡>0; for all 𝑡∈𝐑₀₊ and Φ_𝛼𝑗0(𝑡)=∑^𝑁ℓ=0(𝐴^ℓ0𝑡^ℓ/ℓ!)(𝑡^{(𝛼−1)ℓ}ℓ!/(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑗))=𝐼+𝐴₀𝑡+𝑜(𝑡)>0 for some sufficiently small 𝑡∈𝐑₀₊ for any 𝛼∈𝐑₊.
(b) 𝑁≤∞ and 𝐴₀≥0 since 𝑡^{(𝛼−1)ℓ}ℓ!/(𝛼ℓ+𝑗)Γ(𝛼ℓ+𝑗)≥0; forallℓ∈𝐙₀₊, for all 𝑗∈𝑘−1∪{0} for any 𝛼∈𝐑₊. It follows from inspection of (4.2) since 𝑒^𝑇𝑖(𝑒^𝐴0𝑡)𝑒_𝑚≥0 for all (𝑖,𝑚)∈𝑛×𝑛, since 𝐴₀∈𝑀𝐑^𝑛×𝑛. This implies (Φ_{𝛼𝑗0𝑖𝑚}(𝑡))≥0; for all 𝑡∈𝐑₀₊.
(c) 𝑁<∞⇒min_{0≤ℓ≤𝑁}(𝑡^{(𝛼−1)ℓ}/Γ(𝛼ℓ+𝑘))>0; for all 𝑗∈𝑘−1∪{0} for any 𝛼∈𝐑₊ so that 𝑒^𝑇𝑖(𝑒^𝐴0𝑡)𝑒_𝑚≥0, for all (𝑖,𝑚)∈𝑛×𝑛, since 𝐴₀∈𝑀𝐑^𝑛×𝑛, irrespectively of 𝐴₀≥0 or not, what follows from (4.3). This implies (Φ_{𝛼𝑗0𝑖𝑚}(𝑡))≥0; for all 𝑡∈𝐑₀₊.
(d) 𝑁≤∞, 0<𝛼≤𝑘=1. Then, 𝑗=0 so that 𝑡^{(𝛼−1)ℓ}ℓ!Γ(𝛼ℓ+𝑗+1)=𝑡^{(𝛼−1)ℓ}ℓ!Γ(𝛼ℓ+1)=𝑡^{(𝛼−1)ℓ}ℓ!𝛼ℓΓ(𝛼ℓ)=𝑡^{(𝛼−1)ℓ}(ℓ−1)!𝛼Γ(𝛼ℓ)=𝛼𝑡Γ(^{(1−𝛼)ℓ}ℓΓ)(𝛼ℓ)≥𝛼𝑡1^{(1−𝛼)ℓ}∀ℓ∈𝐙₀₊,,∀𝑡∈𝐑₀₊(4.4) since 0<𝛼≤1 implies Γ(𝛼ℓ)=^∞0𝜏^𝛼ℓ−1𝑒^−𝜏𝑑𝜏≤Γ(ℓ)=^∞0𝜏^ℓ−1𝑒^−𝜏𝑑𝜏,∀ℓ∈𝐙₀₊.(4.5) As a result, Φ_𝛼00(𝑡) from (4.2); for all 𝑡∈𝐑₀₊. Also, direct calculations with (3.3)-(3.4) lead to Φ_𝛼(𝑡)∶=𝑡^𝛼−1𝐸_𝛼𝛼𝐴₀𝑡^𝛼=^∞_ℓ=0𝑡^𝛼−1𝐴^ℓ0𝑡^𝛼ℓΓ((ℓ+1)𝛼)=^∞_ℓ=0𝐴^ℓ0𝑡^ℓℓ!𝑡^{(𝛼−1)(ℓ+1)}ℓ!Γ((ℓ+1)𝛼)(4.6) and similar developments to the above ones yield (Φ_𝛼(𝑡))_𝑖𝑚≥0; for all (𝑖,𝑚)∈𝑛×𝑛, for all 𝑡∈𝐑₀₊ under the same conditions as above in the cases (a) to (d) for Φ_𝛼00(𝑡). On the other hand, one gets from (3.2)–(3.4) for the unforced system with point initial conditions at 𝑡=0: 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Φ_𝛼𝑗0(𝑡)𝑥_𝑗0=Φ_𝛼00(𝑡),…,Φ_{𝛼,𝑘−1,0}(𝑡)𝑥^𝑇00,…,𝑥^{𝑇𝑘−1,0}^𝑇(4.7) which leads to 𝑥_𝛼(𝑡)=Φ_𝛼𝑖0(𝑡)𝑥_𝑖0 by taking point initial conditions 𝑥_𝑖0≠0,𝑥_𝑗0=0,(𝑖≠𝑗),  𝑗∈𝑘−1∪{0} so that Φ_𝛼𝑖0(𝑡) is nonsingular for all 𝑡∈𝐑₀₊ since otherwise the solution is not unique for each given set of initial conditions since any trajectory solution subject to some set of initial conditions 𝑥_𝑖0≠0, 𝑥_𝑗0=0, would have infinitely many initial conditions, subject to identical constraint, so that such a trajectory is not unique which is a contradiction. Since this reasoning may be made for any 𝑗∈𝑘−1∪{0}, Φ_𝛼𝑗0(𝑡) is nonsingular for all 𝑗∈𝑘−1∪{0}, all and, in addition, Φ_𝛼𝑗0(𝑡)>0; forall𝑗∈𝑘−1∪{0}, for all 𝑡∈𝐑₀₊ if either 𝐴₀≥0 or if 𝐴₀ is nilpotent or if 0<𝛼≤𝑘=1 or without these restricting condition within some first interval [0,𝑡). The following properties have been proven:
(a) 𝐴₀∈𝑀𝐑^𝑛×𝑛⇔Φ_𝛼00(𝑡)>0; for all 𝑡∈𝐑₀₊,
(b) 𝐴₀∈𝑀𝐑^𝑛×𝑛⇒(Φ_𝛼𝑗0(𝑡)>0∧detΦ_𝛼𝑗0(𝑡)≠0∧Φ_𝛼(𝑡)≥0; for all 𝑗∈𝑘−1∪{0}, 𝑘=[𝛼]+1 if 𝛼∉𝐙₊and 𝑘=𝛼∈𝐙₊, for all 𝑡∈𝐑₀₊).
It remains to prove Φ_𝛼𝑗0(𝑡)>0; for all 𝑡∈𝐑₀₊⇒𝐴₀∈𝑀𝐑^𝑛×𝑛; for all 𝑗∈𝑘−1, some 𝑡∈𝐑₀₊. This is equivalent to its contrapositive logic proposition. Proceed by contradiction by assuming thereexist𝑗∈𝑘−1 such that 𝐴₀∉𝑀𝐑^𝑛×𝑛⇒Φ_𝛼𝑗0(𝑡)<0, some 𝑡∈𝐑₀₊. Note that 𝐴₀∉𝑀𝐑^𝑛×𝑛⇒𝑒^𝑇𝑖𝑒^𝐴0𝑡(𝑡)𝑒_𝑚=𝑒^𝑇𝑖Φ_𝛼00(𝑡)𝑒_𝑚<0, some 𝑡∈𝐑₀₊, some (𝑖,𝑚≠𝑖)∈𝑛. Then, one gets Φ_{𝛼𝑗0𝑖𝑚}(𝑡)=𝑒^𝑇𝑖Φ_𝛼𝑗0(𝑡)𝑒_𝑚≤𝑒^𝑇𝑖𝑒^𝐴0𝑡𝑒_𝑚max_{0≤ℓ≤∞}𝑡^{(𝛼−1)ℓ}ℓ!Γ(𝛼ℓ+𝑘)<0(4.8) which contradicts (−Φ_𝛼𝑗0(𝑡))>0; for all𝑡∈𝐑₀₊⇐𝐴₀∉𝑀𝐑^𝑛×𝑛; for all 𝑗∈𝑘−1, some 𝑡∈𝐑₀₊. Thus, the proof of Properties (i)-(ii) becomes complete since the above proven property (a) extends to any 𝑗∈𝑘−1∪{0} as follows.
(c) 𝐴₀∈𝑀𝐑^𝑛×𝑛⇔Φ_𝛼𝑗0(𝑡)>0; for all 𝑗∈𝑘−1∪{0}, 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, forall𝑡∈𝐑₀₊; for all 𝛼∈𝐑₊so that the unforced solution for any set of nonnegative point initial conditions is nonnegative for all time and, furthermore, 𝜑_𝑖(𝑡)≥0(∀𝑖∈𝑘−1∪{0}); for all 𝑡∈[−ℎ,0],𝑢(𝑡)∈𝐑ⁿ₀₊; for all 𝑡∈𝐑₀₊,𝐴_𝑖≥0(∀𝑖∈𝑝)and𝐵≥0; for all 𝑡∈𝐑₀₊ implies that (3.2) is everywhere nonnegative within its definition domain. The converse is also true as it follows by contradiction arguments. If there is one entry of B or 𝐴_𝑖(𝑠ome𝑖∈𝑝) which is negative, or if 𝐴₀∉𝑀𝐑^𝑛×𝑛, it can always be found a control 𝑢(𝑡)∈𝐑^𝑛0+ of sufficiently large norm along a given time interval such that some component of the solution is negative for some time. It can be also found that some nonnegative initial condition of sufficiently large norm at 𝑡=0 such that some component of the solution is negative at 𝑡=0⁺. Thus, Property (iii) is proven.

The following result is obvious from the proof of Theorem 4.1.

Corollary 4.2. Theorem 4.1(iii) is satisfied also independent of the delays for any given set of delays satisfying the constraint 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝=ℎ<∞.

Proof. It follows directly since Theorem 4.1 is an independent of the delay size type result and, under the delay constraint 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝=ℎ<∞, it has also to be fulfilled for any combination of delays satisfying the stronger constraint 0=ℎ₀<ℎ₁<ℎ₂<⋯<ℎ_𝑝=ℎ<∞.

Corollary 4.3. Any solution (3.8), subject to (3.9), to the Caputo fractional differential system (3.1) under the delay constraint 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝=ℎ<∞ is nonnegatively independent of the delays within a first interval, that is, it satisfies 𝑥_𝛼(𝑡)∈𝐑^𝑛0+; for all 𝑡∈[−ℎ,𝑡)∩𝐑₀₊ for some sufficiently small 𝑡∈𝐑₀₊ for any given absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0} and any given piecewise continuous vector function 𝑢∶𝐑₀₊→𝐑^𝑛0+ with 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, for all 𝑡∈𝐑₀₊; for all 𝛼∈𝐑₊ if and only if 𝐴₀∈𝑀𝐑^𝑛×𝑛, 𝐴_𝑖∈𝐑^𝑛×𝑛0+(∀𝑖∈𝑝), and 𝐵∈𝐑^𝑛×𝑚0+. In addition, 𝑥_𝛼∶[−ℎ,0)∪𝐑₀₊→𝐑^𝑛0+ if, in addition, either 𝐴₀≥0 or if 𝐴₀ is nilpotent or if 0<𝛼≤𝑘=1. Furthermore, Ψ_𝛼𝑗0(𝑡)>0 (with at least n entries being positive), detΨ_𝛼𝑗0(𝑡)>0(∀𝑗∈𝑘−1∪{0}) and Ψ_𝛼(𝑡)≥0; for all𝑡∈𝐑₀₊(if𝐵∈𝐑^+𝑛×𝑚thenΨ_𝛼(𝑡)>0;∀𝑡∈𝐑₀₊).

Proof. The solution (3.8) is identical to the unique solution (3.2) for (3.1) thus it is everywhere nonnegative under the same conditions that those of Theorem 4.1 which have been extended in Corollary 4.2.

Note that the conditions of nonnegativity of the solution of the above theorem also imply the excitability of all the components of the state-trajectory solution; that is its strict positivity for some 𝑡∈𝐑₊ provided that 𝐵≫0 and the control 𝑢∶𝐑₀₊→𝐑^𝑛0+ is admissible (i.e., piecewise continuous) and nonidentically zero since Ψ_𝛼(𝑡)>0 and nonsingular for all 𝑡∈𝐑₊. It is now seen that the positivity conditions for the Riemann-Liouville fractional differential system (3.5) are not guaranteed in general by the above results for any given absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0} and any given piecewise continuous vector function 𝑢∶𝐑₀₊→𝐑^𝑛0+ with 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, for all 𝑡∈𝐑₀₊; for all 𝛼∈𝐑₊. The following two results hold by using Corollary 3.2 and Corollary 3.4.

Theorem 4.4. Any solution (3.7), subject to (3.3)-(3.4), to the Riemann-Liouville fractional differential system (3.5) under the delay constraint 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝=ℎ<∞ is everywhere nonnegative independent of the delays, that is, it satisfies 𝑥_𝛼∶[−ℎ,0)∪𝐑₀₊→𝐑^𝑛0+, for any given absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0} and any given piecewise continuous vector function 𝑢∶𝐑₀₊→𝐑^𝑛0+ with 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, forall𝑡∈𝐑₀₊; for all 𝛼∈𝐑₊ if 𝐴₀∈𝑀𝐑^𝑛×𝑛, 𝐴_𝑖∈𝐑^𝑛×𝑛0+(∀𝑖∈𝑝), (𝐸_𝛼,𝑗+1(𝑡)−(1/𝑗!)𝐼_𝑛)≥0; forall𝑗∈𝑘−1∪{0}, for all 𝑡∈𝐑₀₊ and 𝐵∈𝐑^𝑛×𝑚0+. The conditions 𝐴₀∈𝑀𝐑^𝑛×𝑛, 𝐴_𝑖∈𝐑^𝑛×𝑛0+(∀𝑖∈𝑝) and 𝐵∈𝐑^𝑛×𝑚0+ are also necessary for 𝑥_𝛼∶[−ℎ,0)∪𝐑₀₊→𝐑^𝑛0+ for any nonnegative function of initial conditions and nonnegative controls. The condition (𝐸_𝛼,𝑗+1(𝑡)−(1/𝑗!)𝐼_𝑛)≥0; for all 𝑗∈𝑘−1∪{0}, for all 𝑡∈𝐑₀₊ is removed for initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+ subject to 𝜑_𝑗(0)=𝑥_𝑗0=0.

Proof. The proof follows in a similar way as the sufficiency part of the proof of Theorem 4.1(iii) by inspecting the nonnegative of the solution Corollary 3.2, (3.7) for a nonnegative function of initial conditions and any nonnegative control.

Theorem 4.5. Any solution (3.10), subject to (3.3)-(3.4), to the Riemann-Liouville fractional differential system (3.5) under the delay constraint 0=ℎ₀≤ℎ₁≤ℎ₂≤⋯≤ℎ_𝑝=ℎ<∞ is everywhere nonnegatively independent of the delays, that is, it satisfies 𝑥_𝛼∶[−ℎ,0)∪𝐑₀₊→𝐑^𝑛0+, for any given absolutely continuous functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+, 𝑗∈𝑘−1∪{0} and any given piecewise continuous vector function 𝑢∶𝐑₀₊→𝐑^𝑛0+ with 𝑘=[𝛼]+1 if 𝛼∉𝐙₊ and 𝑘=𝛼∈𝐙₊, for all 𝑡∈𝐑₀₊; forall𝛼∈𝐑₊ if and only if 𝐴₀∈𝑀𝐑^𝑛×𝑛, 𝐴_𝑖∈𝐑^𝑛×𝑛0+(∀𝑖∈𝑝), (Ψ_𝛼𝑗0(𝑡)−(𝑡^𝑗/𝑗!)𝐼_𝑛)≥0; for all 𝑗∈𝑘−1∪{0}, for all 𝑡∈𝐑₀₊ and 𝐵∈𝐑^𝑛×𝑚0+. The condition (Ψ_𝛼𝑗0(𝑡)−(𝑡^𝑗/𝑗!)𝐼_𝑛)≥0; for all 𝑗∈𝑘−1∪{0}, for all 𝑡∈𝐑₀₊ is removed for initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛0+ subject to 𝜑_𝑗(0)=𝑥_𝑗0=0.

Proof. The proof of sufficiency follows in a similar way as the sufficiency part of the proof of Theorem 4.1(iii) (see also the proof of Theorem 4.5) by inspecting the nonnegativity of the solution Corollary 3.2, (3.7) for a nonnegative function of initial conditions and any nonnegative control. The proof necessity follows by contradiction by inspecting the solution (3.10) as follows.
(a) Assume that 𝐴₀∉𝑀𝐑^𝑛×𝑛 and the solution is nonnegative for all time for any nonnegative function of initial conditions and controls. Take initial conditions 𝜑_𝑗(𝑡)=0; for all 𝑡∈[−ℎ,0], for all 𝑗∈𝑘−1; 𝜑₀(𝑡)=0; for all 𝑡∈[−ℎ,0), 𝜑₀(0)=𝑥₀₀(0)≠0and 𝑢≡0 on 𝐑₀₊. Then (3.10) becomes 𝑥_𝛼(𝑡)=Ψ_𝛼00(𝑡)−𝐼_𝑛𝑥₀₀=𝐸_𝛼00(𝑡)−𝐼_𝑛𝑥₀₀for𝑡∈0,ℎ₁(4.9) since Ψ_𝛼00(𝑡)=Φ_𝛼00(𝑡) for 𝑡∈[0,ℎ₁]. Since 𝐴₀∉𝑀𝐑^𝑛×𝑛, there exist 𝑡∈[0,ℎ₁] and (𝑖=𝑖(𝑡),𝑚(𝑡)≠𝑖)∈𝑛 such that (Φ_𝛼00(𝑡))_𝑖𝑚<0. Otherwise, if 𝐴₀∉𝑀𝐑^𝑛×𝑛 and Φ_𝛼00(𝑡)>0; for all 𝑡∈[0,ℎ₁], it would follow from (4.3) that Φ_𝛼00(𝑡)>0; for all 𝑡∈𝐑₀₊ since 𝑒^𝐴0𝑡=𝑒^{𝜒𝐴0ℎ1+𝐴0𝛿}=𝑒^{𝜒𝐴0ℎ1}𝑒^𝐴0𝛿=𝑒^𝐴0ℎ1^𝜒𝑒^𝐴0𝛿>0(4.10) from the semigroup property of (𝑒^𝐴0𝑡,𝑡∈𝐑₀₊) with 𝜒=𝜒(𝑡,ℎ₁)=[𝑡/ℎ₁] and (0,ℎ₁)∋𝛿=𝛿(𝑡,ℎ₁)=𝑡−𝜒ℎ₁ what implies Φ_𝛼00(𝑡)>0; for all 𝑡∈𝐑₀₊ from (4.3). Thus, 𝐴₀∈𝑀𝐑^𝑛×𝑛 which contradicts 𝐴₀∉𝑀𝐑^𝑛×𝑛. It has been proven that 𝐴₀∉𝑀𝐑^𝑛×𝑛⇒𝑒^𝑇𝑖Φ_𝛼00(𝑡)𝑒_𝑚<0; forall𝑡∈(0,ℎ₁] for some (𝑖=𝑖(𝑡),𝑚(𝑡)≠𝑖)∈𝑛. Now, take 𝑥_00𝑗=𝛿_𝑗𝑚(∀𝑗∈𝑛) where 𝛿_𝑗𝑚 denotes the Kronecker delta. Then, 𝑥_𝛼𝑖(𝑡)=𝑒^𝑇𝑖Ψ^𝛼00(𝑡)−𝑡^𝑗𝑗!𝐼^𝑛𝑥⁰⁰=𝑒^𝑇𝑖Φ^𝛼00(𝑡)−𝑡^𝑗𝑗!𝐼^𝑛𝑒^𝑚=𝑒^𝑇𝑖Φ_𝛼00(𝑡)𝑒_𝑚<0.(4.11) As a result, 𝐴₀∈𝑀𝐑^𝑛×𝑛 is a necessary condition for the solution to be nonnegative for all time irrespective of the delay sizes.
(b) Assume that the solution is nonnegative for all time for any nonnegative function of initial conditions and controls. Assume that 𝑒^𝑇𝑖𝐴_ℓ𝑒_𝑗<0 and ℎ_ℓ≠ℎ_𝑖; for all 𝑖(≠ℓ)∈𝑝 for some 𝑖,𝑗∈𝑛, ℓ∈𝑝. Take initial conditions 𝑥_𝑗0=𝜑_𝑗(0)=0; forall𝑡∈[−ℎ,0]; forall𝑗∈𝑘−1∪{0}, 𝜑_𝑗≡0; forall𝑗∈𝑘−1 and 𝑢≡0. One gets from (3.2) 𝑥_𝛼𝑖(𝑡)=^ℎ′₀𝑒^𝑇𝑖Φ_𝛼(𝑡−𝜏)⎛⎜⎜⎝_{𝑗(≠ℓ)∈𝑝}𝐴_𝑗⎞⎟⎟⎠𝜑₀+𝜏−ℎ^𝑑𝜏_ℎℓ0𝑒^𝑇𝑖Φ_𝛼(𝑡−𝜏)𝐴_ℓ𝜑₀𝜏−ℎ_ℓ𝑑𝜏;∀𝑡∈0,ℎ_𝑖(4.12) for the case ℎ^=ℎ_𝑖; forall𝑖(≠ℓ∈𝑝). Now, if ℎ=ℎ_ℓ>ℎ^, take a further specification of initial conditions as follows: 𝜑₀(𝑡)=0; forall𝑡∈[0,ℎ^], and 𝜑₀(𝜏)=(𝑘₁,…,𝑘_𝑛)^𝑇≫0; forall𝑡∈(ℎ^,ℎ_ℓ] then 𝑥_𝛼𝑖(𝑡)==^ℎℓ_ℎ′𝑒^𝑇𝑖Φ_𝛼(𝑡−𝜏)𝐴_ℓ𝜑₀𝜏−ℎ_ℓ𝑑𝜏_𝑛_𝑛𝑟=1_𝑚=1^ℎℓ_ℎ′Φ_𝛼𝑖𝑟(𝑡−𝜏)𝐴_ℓ𝑟𝑚𝑑𝜏𝑘_𝑚=^𝑛_𝑗=1𝑘_𝑗^ℎℓ_ℎ′Φ^𝑇𝛼𝑖(𝑡−𝜏)𝑑𝜏𝐴_ℓ𝑗=_{(𝑚≠𝑗)∈𝑛}𝑘_𝑚^𝑚_𝑟=1^ℎℓ_ℎ′Φ_𝛼𝑖𝑟(𝑡−𝜏)𝑑𝜏𝐴_ℓ𝑟𝑚+𝑘_𝑗^𝑚_𝑟=1^ℎℓ_ℎ′Φ_𝛼𝑖𝑟(𝑡−𝜏)𝑑𝜏𝐴_ℓ𝑟𝑗,∀𝑡∈0,ℎ_𝑖.(4.13) As a result, 𝐴_𝑖≥0(∀𝑖∈𝑝) is a necessary condition for the solution to be nonnegative for all time irrespective of the delay sizes.
(c) Assume that the solution is nonnegative for all time for any nonnegative function of initial conditions and controls, and 𝐵≥0 is not fulfilled so that it exists at least an entry 𝐵_ℓ𝑗<0 of 𝐵. Then, one has under identically zero initial conditions the following unique solution: 𝑥_𝛼𝑖(𝑡)==^𝑡0𝑒^𝑇𝑖Ψ_𝛼(𝑡−𝜏)𝐵𝑢(𝜏)𝑑𝜏_𝑚_𝑖=1^𝑡0=Ψ^𝑇𝛼𝑖(𝑡−𝜏)𝐵_𝑖𝑢_𝑖(𝜏)𝑑𝜏_𝑚_𝑛𝑖=1_ℓ=1^𝑡0=Ψ_𝛼𝑖ℓ(𝑡−𝜏)𝐵_ℓ𝑖𝑢_𝑖(𝜏)𝑑𝜏_{(𝑖≠𝑗)∈𝑚𝑛}_ℓ=1^𝑡0−Ψ_𝛼𝑖ℓ(𝑡−𝜏)𝐵_ℓ𝑖𝑢_𝑖(𝜏)𝑑𝜏_𝑛_ℓ=1^𝑡0Ψ_𝛼𝑗ℓ(𝑡−𝜏)||𝐵_ℓ𝑗||𝑑𝜏𝑘_𝑢𝑗<0,𝑡∈𝐑₀₊(4.14) provided that 𝑘_𝑢𝑗>∑_{(𝑖≠𝑗)∈𝑚}∑^𝑛ℓ=1∫^𝑡0Ψ_𝛼𝑖ℓ(𝑡−𝜏)𝐵_ℓ𝑖𝑢_𝑖(𝜏)𝑑𝜏/(∑^𝑛ℓ=1∫^𝑡0Ψ_𝛼𝑗ℓ(𝑡−𝜏)|𝐵_ℓ𝑗|𝑑𝜏) by assuming that 𝐵≥0 fails because 𝐵_ℓ𝑗<0 for some (ℓ,𝑗)∈𝑛×𝑚 and a constant control component 𝑢_𝑗≡𝑘_𝑢𝑗>0 is injected on the time interval [0,𝑡] for some arbitrary 𝑡∈𝐑₊ for the remaining control components being chosen 𝑡  be nonnegative for all time. This contradicts that the solution is nonnegative for all time if the condition 𝐵≥0 fails.

Remark 4.6. Note that Theorem 4.1 can be extended as a necessary condition for 𝑡∈[0,ℎ₁] since Ψ_𝛼𝑗0(𝑡)=Φ_𝛼𝑗0(𝑡) for 𝑡∈[0,ℎ₁]; forall𝑗∈𝑘−1∪{0}, forall𝑡∈𝐑₀₊.

Remark 4.7. Note by simple calculation that (𝐴₀∈𝑀𝐑^𝑛×𝑛∧𝐴_𝑖≥0(𝑖∈𝑝))⇒(∑^𝑝𝑖=0𝐴_𝑖)∈𝑀𝐑^𝑛×𝑛. This is a necessary and sufficient condition for the nonnegativity of the solutions of the Caputo fractional differential system (3.1) of arbitrary order 𝛼∈𝐑₊ under arbitrary nonnegative controls and initial conditions in the absence of delays; that is, for ℎ_𝑖=0; forall𝑖∈𝜔∪{0} and any 𝜔∈𝐙₊.

Remark 4.8. The given conditions to guarantee that the solution is everywhere nonnegative under any given arbitrary nonnegative initial conditions and nonnegative controls are independent of the sizes of the delays type; that is, for any given set of 𝑝 delays. However, the conditions are weakened for particular situations involving repeated delays as follows. Note from Theorem 4.5 that the various given conditions 𝐴_𝑖≥0 of necessary type to guarantee the nonnegativity of the solution under any admissible nonnegative controls and nonnegative initial conditions are weakened to (∑^𝜈𝑖_ℓ=1𝐴_(∑^{𝑖−1𝑗=0}_{𝜈𝑗+ℓ)}) if there is some repeated delay ℎ_𝑖 of multiplicity 𝜈_𝑖≥2 (i.e., the number of distinct delays is 0≤𝑞<𝑝). Also, if ℎ₀=0 is repeated with multiplicity 𝜈₀≥2 then the condition 𝐴₀∈𝑀𝐑^𝑛×𝑛 for 𝜈₀=1 is replaced by (∑^𝜈0−1_ℓ=0𝐴_ℓ)∈𝑀𝐑^𝑛×𝑛.

Remark 4.9. Note that there is a duality of all the given results of sufficiency type or necessary and sufficiency type in the sense that the solutions are guaranteed to be nonpositive for all time under similar conditions for the cases when all components of the controls and initial conditions are nonpositive for all time.

5. Asymptotic Behavior of Unforced Solutions for 𝛼∈𝐑₊

The asymptotic behaviour and the stability properties of the Caputo fractional differential system (3.1) can be investigated via the extension of the subsequent formulas for 𝛼∈𝐑₊, (see (1.8.27)–(1.8.29), [1]).(1)If 0<𝛼<2 then for |𝑧|→∞ and some 𝜇∈𝐑 satisfying 𝜇<𝜋min(1,𝛼):𝐸_𝛼𝛽(𝑧)=1𝛼𝑧^{(1−𝛽)/𝛼}𝑒^(𝑧1/𝛼)−^𝑁_𝑗=11Γ(𝛽−𝛼𝑗)𝑧1^𝑗+𝑂1𝑧^𝑁+1(5.1) with |arg𝑧|≤𝜇<𝜋min(1,𝛼), any 𝑁∈𝐙₊, and 𝐸_𝛼𝛽(𝑧)=−^𝑁_𝑗=11Γ(𝛽−𝛼𝑗)𝑧1^𝑗+𝑂1𝑧^𝑁+1(5.2) with 𝜋≥|arg𝑧|≥𝜇<𝜋=𝜋min(1,𝛼), any 𝑁∈𝐙₊.(2)If 𝛼≥2 then for |𝑧|→∞𝐸_𝛼𝛽(𝑧)=1𝛼_𝑗∈𝑄𝑧^1/𝛼𝑒^{2𝑗𝜋𝐢/𝛼}^1−𝛽𝑒^{(𝑒2𝑗𝜋𝐢/𝛼𝑧1/𝛼)}−^𝑁_𝑗=11Γ(𝛽−𝛼𝑗)𝑧1^𝑗+𝑂1𝑧^𝑁+1(5.3) for any 𝑁∈𝐙₊ with 𝛽∈𝑘−1∪{0}∪{𝛼}, |arg𝑧|≤𝛼𝜋/2, 𝑄∶={𝑛∈𝐙∶|arg𝑧+2𝜋𝑛|≤𝛼𝜋/2} and 𝐢=√−1 being the complex imaginary unit. The above formulas are extendable to the Mittag-Leffler matrix functions 𝐸_𝛼,𝑗+1(𝐴₀𝑡^𝛼)∶=∑^∞ℓ=0(𝐴₀𝑡^𝛼)^ℓ/Γ(𝛼ℓ+𝑗+1); forall𝑗∈𝑘−1∪{0}, respectively, 𝐸_𝛼𝛼(𝐴₀𝑡^𝛼) by identifying 𝑧→𝐴₀𝑡^𝛼, 𝑧^1/𝛼→(𝐴₀)^1/𝛼𝑡 (if (𝐴₀)^1/𝛼exists) and 𝑧⁻¹→𝐴⁰⁻¹𝑡^−𝛼 (if 𝐴₀ is non-singular), 𝛽→𝑗+1, respectively, 𝛽→𝛼. Irrespective of the existence of (𝐴₀)^1/𝛼and of 𝐴₀≠0 being singular or nonsingular, it is possible to identify 𝑧⁻¹→‖𝐴₀‖⁻¹𝑡^−𝛼 and 𝑧→‖𝐴₀‖𝑡^𝛼 and to use‖‖𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼‖‖≤𝐸_𝛼,𝑗+1‖‖𝐴₀‖‖𝑡^𝛼=^∞_ℓ=0‖‖𝐴₀‖‖𝑡^𝛼^ℓΓ(𝛼ℓ+𝑗+1),∀𝑗∈‖‖𝐸𝑘−1∪{0},_𝛼𝛼𝐴₀𝑡^𝛼‖‖≤𝐸_𝛼𝛼‖‖𝐴₀‖‖𝑡^𝛼=^∞_ℓ=0‖‖𝐴₀‖‖𝑡^𝛼^ℓΓ(𝛼ℓ+𝛼).(5.4)

The method may be used to calculate an asymptotic estimate of the solution (3.2) if 𝐴₀ is non-singular (or an upperbounding function for any nonzero 𝐴₀) of the Caputo fractional differential system (3.1), via (3.3)-(3.4), or, equivalently (3.8), via (3.9) and (3.3)-(3.4). The estimations may be extended with minor modification to the Riemann-Liouville fractional differential system (3.5). Note that if all the complex eigenvalues of 𝐴₀ appear by conjugate pairs 𝐴₀ then 𝐴₀=𝑇⁻¹𝐽₀𝐴₀ where 𝐽₀ is its real canonical form. First, consider two separate cases as follows.

(A) Assume that 𝛼∈𝐑₊, 𝐴₀ is real non-singular and (𝐴₀)^1/𝛼 exists; that is, there exist 𝑀 such that 𝑀^𝛼=𝐴₀ and 𝐴_𝑖(𝑖∈𝑝) is real. Then, one gets from (5.1)–(5.3):𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼=1𝛼𝐴^0−𝑗^1/𝛼𝑡^−𝑗𝑒^{(𝐴01/𝛼)𝑡}−^𝑁_ℓ=11Γ((1−𝛼)ℓ+1)𝐴^0−ℓ𝑡^−ℓ𝛼𝐸+𝑂𝐴^0−(𝑁+1)𝑡^{−(𝑁+1)𝛼},∀𝑗∈𝑘−1∪{0},_𝛼𝛼𝐴₀𝑡^𝛼=𝛼1𝐴^{0(1−𝛼)/𝛼}𝑡^1−𝛼𝑒^{(𝐴01/𝛼)𝑡}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝑡^−ℓ𝛼+𝑂𝐴^0−(𝑁+1)𝑡^{−(𝑁+1)𝛼}(5.5) as 𝑡→∞ if 0<𝛼<2, for any 𝑁∈𝐙₊, 𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼=1𝛼_ℓ∈𝑄𝐴^01/𝛼𝑡𝑒^{2ℓ𝜋𝐢/𝛼}^−𝑗𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝑡)}−^𝑁_ℓ=11Γ((1−𝛼)ℓ+1)𝐴^0−ℓ𝑡^−ℓ𝛼𝐸+𝑂𝐴^0−(𝑁+1)𝑡^{−(𝑁+1)𝛼},∀𝑗∈𝑘−1∪{0},_𝛼𝛼𝐴₀𝑡^𝛼=𝛼1_ℓ∈𝑄𝐴^01/𝛼𝑡𝑒^{2ℓ𝜋𝐢/𝛼}^1−𝛼𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝑡)}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝑡^−ℓ𝛼+𝑂𝐴^0−(𝑁+1)𝑡^{−(𝑁+1)𝛼}(5.6) as 𝑡→∞ if 𝛼≥2, for any 𝑁∈𝐙₊, with 𝑄∶={𝑛∈𝐙∶|𝑛|≤𝛼/4}.

(B) Assume that 𝛼∈𝐑₊ and 𝐴_𝑖(𝑖∈𝑝∪{0}) is real, one obtains from (5.1)-(5.2):||𝐸_𝛼𝛽(𝑧)||≤𝐸_𝛼𝛽(𝑧)=1𝛼|𝑧|^{(1−𝛽)/𝛼}||𝑒^𝑧||^1/𝛼+|||||^𝑁_ℓ=11Γ(𝛽−𝛼ℓ)𝑧1^𝑗+𝑂1𝑧^𝑁+1|||||(5.7) for 0<𝛼<2, 𝛽∈𝑘−1∪{0},||𝐸_𝛼𝛽(𝑧)||≤𝐸_𝛼𝛽(𝑧)=1𝛼_𝑗∈𝑄|𝑧|^{(1−𝛽)/𝛼}||𝑒^𝑧||^1/𝛼+|||||^𝑁_ℓ=11Γ(𝛽−𝛼ℓ)𝑧1^𝑗+𝑂1𝑧^𝑁+1|||||(5.8) for 𝛼≥2, 𝛽∈𝑘−1∪{0,𝛼}. Thus, on gets from (5.7)‖‖𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼‖‖≤1𝛼‖‖𝐴⁰‖‖^−𝑗^1/𝛼𝑡^−𝑗‖‖𝑒^𝐴0𝑡‖‖^1/𝛼−^𝑁_ℓ=11Γ((1−𝛼)ℓ+1)‖‖𝐴⁰‖‖^−ℓ𝑡^−ℓ𝛼‖‖𝐸+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{−(𝑁+1)𝛼},∀𝑗∈𝑘−1∪{0},_𝛼𝛼𝐴₀𝑡^𝛼‖‖≤𝛼1‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}𝑡^1−𝛼‖‖𝑒^𝐴0𝑡‖‖^1/𝛼−^𝑁_ℓ=11Γ((1−ℓ)𝛼)‖‖𝐴⁰‖‖^−ℓ𝑡^−ℓ𝛼+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{−(𝑁+1)𝛼}(5.9) as 𝑡→∞, for any 𝑁∈𝐙₊, if 0<𝛼<2, and one gets from (5.8)‖‖𝐸_𝛼,𝑗+1𝐴₀𝑡^𝛼‖‖≤1𝛼_ℓ∈𝑄‖‖𝐴⁰‖‖^1/𝛼𝑡^−𝑗‖‖𝑒^𝐴0𝑡‖‖^1/𝛼−^𝑁_ℓ=11Γ((1−𝛼)ℓ+1)‖‖𝐴⁰‖‖^−ℓ𝑡^−ℓ𝛼‖‖𝐸+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{−(𝑁+1)𝛼},∀𝑗∈𝑘−1∪{0},_𝛼𝛼𝐴₀𝑡^𝛼‖‖≤𝛼1_ℓ∈𝑄‖‖𝐴⁰‖‖^1/𝛼𝑡^1−𝛼‖‖𝑒^𝐴0𝑡‖‖^1/𝛼−^𝑁_ℓ=11Γ((1−ℓ)𝛼)‖‖𝐴⁰‖‖^−ℓ𝑡^−ℓ𝛼+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{−(𝑁+1)𝛼}(5.10) as 𝑡→∞ if 𝛼≥2, for any 𝑁∈𝐙₊, with 𝑄∶={𝑛∈𝐙∶|𝑛|≤𝛼/4}. The formula (3.8) for the solution is more useful than its equivalent expression (3.2) to investigate the asymptotic properties of the Caputo fractional differential system. Therefore, we obtain now either explicit or upperbounding asymptotic expressions for (3.9) by using (5.5) to (5.9) as follows.

(1) Assume that 𝛼∈𝐑₊, 𝐴₀ is real non-singular, (𝐴₀)^1/𝛼 exists and 𝐴_𝑖(𝑖∈𝑝) are also real. Then, one gets from (5.5)–(5.6) into (3.9):Ψ_𝛼𝑗0(𝑡)=1𝛼𝐴^0−𝑗^1/𝛼𝑒^{(𝐴01/𝛼)𝑡}−^𝑁_ℓ=11Γ(𝑗+1−𝛼ℓ)𝐴^0−ℓ𝑡^{𝑗−ℓ𝛼}++𝑂𝐴^0−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}(5.11)_𝑝_𝑖=1^𝑡0𝛼1𝐴^{0(1−𝛼)/𝛼}𝑒^{(𝐴01/𝛼)𝜏}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝜏^{(1−ℓ)𝛼−1}×𝐴_𝑖+𝑂𝐴^0−(𝑁+1)𝜏^{−𝑁𝛼−1}ΨΨ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖𝑑𝜏,_𝛼(𝑡)=1𝛼𝐴^{0(1−𝛼)/𝛼}𝑒^{(𝐴01/𝛼)𝑡}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝑡^{(1−ℓ)𝛼}+_𝑝+𝑂𝐴^0−(𝑁+1)𝑡^−𝑁𝛼_𝑖=1^𝑡0𝛼1𝐴^{0(1−𝛼)/𝛼}𝑒^{(𝐴01/𝛼)𝜏}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝜏^{(1−ℓ)𝛼−1}×𝐴_𝑖+𝑂𝐴^0−(𝑁+1)𝜏^{−𝑁𝛼−1}Ψ_𝛼𝑡−𝜏−ℎ_𝑖𝑑𝜏(5.12) for all 𝑗∈𝑘−1∪{0} as 𝑡→∞ if 0<𝛼<2, for any 𝑁∈𝐙₊, andΨ_𝛼𝑗0(𝑡)=1𝛼_ℓ∈𝑄𝐴^01/𝛼𝑒^{2ℓ𝜋𝐢/𝛼}^−𝑗𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝑡)}−^𝑁_ℓ=11Γ(𝑗+1−𝛼ℓ)𝐴^0−ℓ𝑡^{𝑗−ℓ𝛼}++𝑂_𝑝𝐴^0−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}_𝑖=1^𝑡0⎛⎜⎜⎜⎝1𝛼_ℓ∈𝑄𝐴^{0(1−𝛼)/𝛼}𝑒^{2ℓ𝜋𝐢/𝛼}^−𝑗𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝜏)}−^𝑁_ℓ=11Γ((1−ℓ)𝛼)𝐴^0−ℓ𝜏^{(1−ℓ)𝛼−1}×𝐴_𝑖+𝑂𝐴^0−(𝑁+1)𝜏^{−𝑁𝛼−1}⎞⎟⎟⎟⎠ΨΨ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖𝑑𝜏,_𝛼(𝑡)=1𝛼_ℓ∈𝑄𝐴^{0(1−𝛼)/𝛼}𝑒^{2ℓ𝜋(1−𝛼)𝐢/𝛼}𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝑡)}−^𝑁_ℓ=11Γ((1−𝛼)ℓ+1)𝐴^0−ℓ𝑡^{𝑗−ℓ𝛼}+_𝑝+𝑂𝐴^0−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}_𝑖=1^𝑡0⎛⎜⎜⎜⎝1𝛼_ℓ∈𝑄𝐴^{0(1−𝛼)/𝛼}𝑒^{2ℓ𝜋(1−𝛼)𝐢/𝛼}𝑒^{(𝑒2ℓ𝜋𝐢/𝛼(𝐴01/𝛼)𝜏)}−^𝑁_ℓ=11Γ((1−𝛼)ℓ)𝐴^0−ℓ𝜏^{(1−ℓ)𝛼−1}×𝐴_𝑖+𝑂𝐴^0−(𝑁+1)𝜏^{−𝑁𝛼−1}⎞⎟⎟⎟⎠Ψ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖𝑑𝜏(5.13) for all 𝑗∈𝑘−1∪{0} as 𝑡→∞ if 𝛼≥2, for any 𝑁∈𝐙₊.

(2) Assume that 𝛼∈𝐑₊ and 𝐴_𝑖(𝑖∈𝑝∪{0}) are real. Then, ‖‖Ψ_𝛼𝑗0(𝑡)‖‖≤1𝛼‖‖𝐴⁰‖‖^−𝑗^1/𝛼‖‖𝑒^𝐴0𝑡‖‖^1/𝛼+^𝑁_ℓ=11||Γ(𝑗+1−𝛼ℓ)||‖‖𝐴⁰+_𝑝‖‖^−ℓ𝑡^{𝑗−ℓ𝛼}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}_𝑖=1^𝑡0⎛⎜⎜⎝1𝛼‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}‖‖𝑒^𝐴0𝜏‖‖^1/𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰‖‖^−ℓ𝜏^{(1−ℓ)𝛼−1}‖‖Ψ+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝜏^{−𝑁𝛼−1}⎞⎟⎟⎠‖‖𝐴^𝑖‖‖‖‖Ψ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖‖‖𝑑𝜏,(5.14)_𝛼(𝑡)‖‖≤1𝛼‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}‖‖𝑒^𝐴0𝑡‖‖^1/𝛼𝑡^𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰+_𝑝‖‖^−ℓ𝑡^{(1−ℓ)𝛼}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^−𝑁𝛼_𝑖=1^𝑡0𝛼1‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}𝜏^𝛼‖‖𝑒^𝐴0𝜏‖‖^1/𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰‖‖^−ℓ𝜏^{(1−ℓ)𝛼−1}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝜏^{−𝑁𝛼−1}‖‖𝐴^𝑖‖‖‖‖Ψ_𝛼𝑡−𝜏−ℎ_𝑖‖‖𝑑𝜏(5.15) for all 𝑗∈𝑘−1∪{0} as 𝑡→∞ if 0<𝛼<2, for any 𝑁∈𝐙₊, and‖‖Ψ_𝛼𝑗0(𝑡)‖‖≤1𝛼‖‖𝐴⁰‖‖^1/𝛼^−𝑗‖‖𝑒^𝐴0𝑡‖‖^1/𝛼+^𝑁_ℓ=11||Γ(𝑗+1−𝛼ℓ)||‖‖𝐴⁰+_𝑝‖‖^−ℓ𝑡^{𝑗−ℓ𝛼}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}_𝑖=1^𝑡0𝛼1‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}‖‖𝑒^𝐴0𝜏‖‖^1/𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰‖‖^−ℓ𝜏^{(1−ℓ)𝛼−1}‖‖Ψ+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝜏^{−𝑁𝛼−1}‖‖𝐴^𝑖‖‖‖‖Ψ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖‖‖𝑑𝜏,_𝛼(𝑡)‖‖≤1𝛼‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}‖‖𝑒^𝐴0𝑡‖‖^1/𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰+_𝑝‖‖^−ℓ𝑡^{𝑗−ℓ𝛼}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝑡^{𝑗−(𝑁+1)𝛼}_𝑖=1^𝑡0𝛼1‖‖𝐴⁰‖‖^{(1−𝛼)/𝛼}‖‖𝑒^𝐴0𝜏‖‖^1/𝛼+^𝑁_ℓ=11||Γ((1−ℓ)𝛼)||‖‖𝐴⁰‖‖^−ℓ𝜏^{(1−ℓ)𝛼−1}+𝑂‖‖𝐴₀‖‖^−(𝑁+1)𝜏^{−𝑁𝛼−1}‖‖𝐴^𝑖‖‖‖‖Ψ_𝛼𝑗0𝑡−𝜏−ℎ_𝑖‖‖𝑑𝜏(5.16) for all 𝑗∈𝑘−1∪{0} as 𝑡→∞ if 𝛼≥2, for any 𝑁∈𝐙₊.

For further discussion, note that there exists a set of linearly independent continuously differential real functions {𝛼_𝑖∶𝐑₀₊→𝐑,𝑖∈𝜈−1∪{0}}, where 𝜈 is the degree of the minimal polynomial of any square real matrix 𝐴₀ such that:𝑒^𝐴0𝑡=^𝜈−1_𝑖=0𝛼_𝑖(𝑡)𝐴^𝑖0=^𝜈_𝜈𝑗=0^𝑖−1_𝑖=0𝑘_𝑖𝑗𝑡^𝑖𝑒^𝜆𝑗𝑡;∀𝑡∈𝐑₀₊(5.17) (see, e.g., [4, 5]), where 𝑘_𝑖𝑗∈𝐑; 𝑖∈𝜈−1∪{0}, 𝑗∈𝜈∪{0}, 𝜎(𝑀)∶={𝜆_𝑖∈𝐂∶det(𝜆_𝑖𝐼_𝑛−𝐴₀)=0} is the spectrum of 𝐴₀ defined by the set of eigenvalues 𝜆_𝑖of 𝑀 of respective index 𝜈_𝑖 (i.e., the multiplicity of 𝜆_𝑖in the minimal polynomial of 𝐴₀) and algebraic multiplicity 𝜇_𝑖 (i.e., the multiplicity of 𝜆_𝑖 in the characteristic polynomial of 𝐴₀) so that 𝑛=∑^𝑛𝑖=1𝑛_𝑖≥𝜈=∑^𝑛𝑖=1𝜈_𝑖 with 𝑛 being the order of 𝐴₀ with 𝜈 being the degree of its minimal polynomial. The subsequent fractional calculus-related stability result is based on the above formulas.

Theorem 5.1. The following properties hold.
(i) If 𝑘=𝛼=1 (the particular standard bon-fractional case) then (3.1) is globally Lyapunov stable independent of the delays if ‖‖‖1𝛽₁𝐴₁,1𝛽₂𝐴₂,…,1𝛽_𝑝𝐴_𝑝‖‖‖₂≤−𝜇₂𝐴₀(5.18) requiring for the ℓ₂-matrix measure of 𝐴₀ to fulfil 𝜇₂(𝐴₀)∶=(1/2)𝜆_max(𝐴₀+𝐴^𝑇0)≤0, for some 𝛽_𝑖∈𝐑₊(𝑖∈𝑝) subject to ∑^𝑝𝑖=1𝛽^2𝑖=1, [6]. Also, Ψ₁₀₀(𝑡)=+𝑒^𝐴0𝑡_𝑝−𝐴⁰⁻¹𝑡⁻¹+𝑂𝐴⁰⁻¹𝑡⁻¹_𝑖=1^𝑡0𝑒^𝐴0𝜏−𝐴⁰⁻¹𝜏⁻¹+𝑂𝐴⁰⁻¹𝜏⁻¹𝐴_𝑖Ψ₁₀₀𝑡−𝜏−ℎ_𝑖𝑑𝜏as𝑡⟶∞(5.19) is bounded provided that 𝐴₀ is non-singular with 𝑒^𝐴0𝑡 being of the form (5.17) if (5.18) holds and then the unforced solution: 𝑥_𝛼(𝑡)=^𝑘−1_𝑗=0Ψ^𝛼𝑗0(𝑡)𝑥_𝑗0+^𝑝_𝑖=1^ℎ𝑖₀Ψ_𝛼𝑗0(𝑡−𝜏)𝜑_𝑗𝜏−ℎ_𝑖𝑑𝜏.(5.20) Is bounded for all time. Furthermore, ‖‖Ψ₁₀₀(𝑡)‖‖≤+‖‖𝑒^𝐴0𝑡_𝑝‖‖+‖‖𝐴₀‖‖⁻¹𝑡⁻¹+𝑂‖‖𝐴₀‖‖⁻²𝑡⁻²_𝑖=1^𝑡0‖‖𝑒^𝐴0𝜏‖‖+‖‖𝐴₀‖‖⁻¹𝜏⁻¹+𝑂‖‖𝐴₀‖‖⁻¹𝜏⁻¹×‖‖𝐴_𝑖‖‖‖‖Ψ₁₀₀𝑡−𝜏−ℎ_𝑖‖‖𝑑𝜏as𝑡⟶∞(5.21) if (5.18) holds irrespective of 𝐴₀ being singular or non-singular. If, in addition, 𝜇₂(𝐴₀)<0 and (5.18) holds with strict inequality then (3.1) is globally asymptotically Lyapunov stable independent of the delays and Ψ₁₀₀(𝑡)⟶^𝑝_𝑖=1^𝑡0𝑒^𝐴0𝜏−𝐴⁰⁻¹𝜏⁻¹+𝑂𝐴⁰⁻¹𝜏⁻¹𝐴_𝑖Ψ₁₀𝑡−𝜏−ℎ_𝑖𝑑𝜏⟶0as𝑡⟶∞.(5.22)
(ii) If 𝑘=1 and 𝛼∈(0,1] the inequality (5.18) is strict then (3.1) is globally Lyapunov stable independent of the delays if 𝜇₂(𝐴^01/𝛼)≤0 and ‖‖‖1𝛽₁𝐴₁,1𝛽₂𝐴₂,…,1𝛽_𝑝𝐴_𝑝‖‖‖₂≤−𝜇₂𝐴^01/𝛼(5.23) provided that 𝐴₀ is non-singular and 𝐴^01/𝛼 exists. Also, then (3.1) is globally asymptotically Lyapunov stable independent of the delays if, in addition, 𝜇₂(𝐴^01/𝛼)<0 and ‖‖‖1𝛽₁𝐴₁,1𝛽₂𝐴₂,…,1𝛽_𝑝𝐴_𝑝‖‖‖₂Ψ<||𝜇₂𝐴^01/𝛼||,(5.24)_𝛼00(𝑡)=1𝛼𝑒^{(𝐴01/𝛼)𝑡}−1Γ(1−𝛼)𝐴⁰⁻¹+_𝑝𝑡^−𝛼+𝑂𝐴⁰⁻²𝑡^−2𝛼_𝑖=1^𝑡0𝛼1𝐴^{0(1−𝛼)/𝛼}×𝐴_𝑖𝑒^{(𝐴01/𝛼)𝜏}−𝐴⁰⁻¹𝜏⁻¹+𝑂𝐴^0−(𝑁+1)𝜏^−𝛼−1Ψ_𝛼00𝑡−𝜏−ℎ_𝑖𝑑𝜏⟶0as𝑡⟶∞.(5.25) If either 𝐴₀ is singular or 𝐴^01/𝛼 does not exists then (5.25) is replaced by a corresponding less than or equal to relation of norms with the replacements 𝐴₀→‖𝐴₀‖, 𝐴⁰⁻¹→‖𝐴₀‖⁻¹, and 𝑒^𝐴0𝑡→‖𝑒^𝐴0𝑡‖.
(iii) Assume that 𝐽_𝐴0=𝐽_𝐴0𝑑+𝐽_𝐴0 is the canonical real form of 𝐴₀ (in particular, its Jordan form if all the eigenvalues are real) with 𝐽_𝐴0𝑑 being diagonal and 𝐽_𝐴0being off diagonal such that the above decomposition is unique with 𝐴₀=𝑇⁻¹𝐽_𝐴0𝑇 where 𝑇 is a unique non-singular transformation matrix. Then, the Caputo fractional differential system (3.1) is globally Lyapunov stable independently of 𝐴^01/𝛼 to exist or not by replacing 𝜇₂(𝐴^01/𝛼)→𝜇₂(𝐽^𝐴1/𝛼_0𝑑) in (5.23) by ‖‖‖1𝛽₀𝑇⁻¹𝐽_𝐴0𝑇,1𝛽₁𝑇⁻¹𝐴₁𝑇,1𝛽₂𝑇⁻¹𝐴₂𝑇,…,1𝛽_𝑝𝑇⁻¹𝐴_𝑝𝑇‖‖‖₂≤||𝜇₂𝐽_𝐴0𝑑||^1/𝛼(5.26) with 𝜇₂(𝐽^𝐴1/𝛼_0𝑑)≤0 for some set of numbers 𝛽_𝑖∈𝐑₊(𝑖∈𝑝∪{0}) satisfying ∑^𝑝𝑖=0𝛽^2𝑖=1. The fractional system is globally asymptotically Lyapunov stable for one such a set of real numbers if 𝜇₂(𝐽^𝐴1/𝛼_0𝑑)<0, what implies that |arg𝜆|<𝛼𝜋/2, for all 𝜆∈𝜎(𝐴₀), and ‖‖‖1𝛽₀𝑇⁻¹𝐽_𝐴0𝑇,1𝛽₁𝑇⁻¹𝐴₁𝑇,1𝛽₂𝑇⁻¹𝐴₂𝑇,…,1𝛽_𝑝𝑇⁻¹𝐴_𝑝𝑇‖‖‖₂<||𝜇₂𝐽_𝐴0𝑑||^1/𝛼.(5.27)

Proof. It turns out that 𝑥_𝛼(𝑡) is bounded for all time so that (3.1) is globally Lyapunov stable if ‖Ψ_𝛼𝑗0(𝑡)‖ is bounded; forall𝑗∈𝑘−1∪{0} for all 𝑡∈𝐑₀₊ for any bounded functions of initial conditions 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛; for all 𝑗∈𝑘−1∪{0} with 𝜑_𝑗(0)=𝑥_𝑗(0)=𝑥_𝑗0. If, in addition, ‖Ψ_𝛼𝑗0(𝑡)‖→0 as 𝑡→∞ then 𝑥_𝛼(𝑡)→0 as 𝑡→∞ so that (3.1) is globally asymptotically Lyapunov stable and the solution (5.20) is bounded for all time. Thus, if 𝑘=𝛼=1 (the particular standard bon-fractional case) then (3.1) is globally Lyapunov stable if ‖Ψ₁₀₀(𝑡)‖ is bounded for all 𝑡∈𝐑₀₊. A sufficient condition independent of the delays is that (5.18) holds requiring trivially for the ℓ₂-matrix measure of 𝐴₀ to fulfil 𝜇₂(𝐴₀)∶=(1/2)𝜆_max(𝐴₀+𝐴^𝑇0)≤0, where the for some 𝛽_𝑖∈𝐑₊(𝑖∈𝑝) subject to ∑^𝑝𝑖=1𝛽^2𝑖=1, [27]. Equation (5.19) follows from (5.11) after inspection for 𝑁=1 and it is bounded as 𝑡→∞ and since otherwise the global stability property (5.18) would fail contradicting its sufficient condition for 𝑗+1=𝑘=𝛼=1. Equation (5.20) follows from (5.14) for 𝑗+1=𝑘=𝛼=𝑁=1 irrespective of 𝐴₀ being singular or non-singular and of the fact of 𝐴^01/𝛼 to exist or not. Equation (5.21) follows from (5.19) since 𝜇₂(𝐴₀)<0 implies that 𝐴₀ is a stability matrix then Re(𝜆)<0; forall𝜆∈𝜎(𝐴₀)  and, furthermore, Ψ₁₀₀(𝑡)→0, and the unforced solution 𝑥_𝛼(𝑡)→0, as 𝑡→∞ from the strict inequality guaranteeing global asymptotic stability independent of the delays, namely, ‖(1/𝛽₁)𝐴₁,(1/𝛽₂)𝐴₂,…,(1/𝛽_𝑝)𝐴_𝑝‖₂<|𝜇₂(𝐴₀)|. Property (i) has been proven. Property (ii) has a similar proof for 𝛼∈(0,1], 𝑘=1 by replacing 𝐴₀→𝐴^01/𝛼. Property (iii) follows by using the matrix similarity transformation 𝐴₀=𝑇⁻¹𝐽_𝐴0𝑇=𝑇⁻¹(𝐽_𝐴0𝑑+𝐽_𝐴0)𝑇 and using the homogeneous transformed Caputo fractional differential system from (3.1): ^𝐶_𝐶𝐷^𝛼0+𝑧(𝑡)=^𝐶𝐷^𝛼0+𝑇𝑥(𝑡)=^𝑝_𝑖=0𝐴_𝑖𝑇𝑥𝑡−ℎ_𝑖⟺𝐷^𝛼0+𝑥(𝑡)=^𝑝_𝑖=0𝑇⁻¹𝐴_𝑖𝑇𝑥𝑡−ℎ_𝑖=𝑇⁻¹𝐴₀𝑇𝑥(𝑡)+^𝑝_𝑖=1=𝑇₋₁𝑇⁻¹𝐴_𝑖𝑇𝑥𝑡−ℎ_𝑖𝐽_𝐴0𝑑𝑇𝑥(𝑡)+^𝑝_𝑖=0𝑇⁻¹𝐴_𝑖𝑇𝑥𝑡−ℎ_𝑖,(5.28) where 𝑧(𝑡)=𝑇𝑥(𝑡); for all 𝑡∈𝐑₀₊, ℎ₀=0 plays the role of an additional delay. 𝐴₀=𝐽_𝐴0 and 𝐴_𝑖=𝐴_𝑖(𝑖∈𝑝) by noting also that since (𝐽_𝐴0𝑑+𝐽^∗𝐴_0𝑑) is diagonal with real eigenvalues by construction, one has |||𝜇₂𝐽^𝐴1/𝛼_0𝑑|||==|||||R12e𝜆𝜆^max_1/𝛼max𝐽^𝐴1/𝛼_0𝑑+𝐽^𝐴1/𝛼_0𝑑^∗|||=|||𝜆^max𝐽^𝐴1/𝛼_0𝑑|||=|||Re𝜆^max𝐽^𝐴1/𝛼_0𝑑|||𝐽_𝐴0𝑑||=||Re𝜆^1/𝛼max𝐴_0𝑑||=||𝜇₂𝐽_𝐴0𝑑||^1/𝛼.(5.29) Then, the proof is similar to that of the related part of Property (ii). Note also that 𝜇₂(𝐽^𝐴1/𝛼_0𝑑)<0, implies that Re𝜆^1/𝛼<0⟺arg𝜆𝛼∈−𝜋2,𝜋2;∀𝜆∈𝜎𝐴⁰⟺||arg𝜆||>𝛼2𝜋,∀𝜆∈𝜎𝐴⁰.(5.30)

Remark 5.2. Note that a similar expressions to (5.25) applies to guarantee global asymptotic stability for 𝛼∈(0,1] in Theorem 5.1(iii) by replacing 𝐴₀→𝑇⁻¹𝐽_𝐴0𝑑𝑇 and 𝐴_𝑖→𝑇⁻¹𝐴_𝑖𝑇 with 𝐴_𝑖(𝑖∈𝑝∪{0}) defined in the proof of Theorem 5.1(iii). Theorem 5.1 establishes that for any stability matrix 𝐴₀, the asymptotic stability condition of sufficient type is as follows:‖‖‖1𝛽₀𝑇⁻¹𝐽_𝐴0𝑇,1𝛽₁𝑇⁻¹𝐴₁𝑇,1𝛽₂𝑇⁻¹𝐴₂𝑇,…,1𝛽_𝑝𝑇⁻¹𝐴_𝑝𝑇‖‖‖₂<||𝜇₂𝐽_𝐴0𝑑||^1/𝛼(5.31) provided that 𝜇₂(𝐽^1/𝛼0_𝐴0𝑑)<0 extends from 𝛼=𝛼₀≤1, (in particular, from the standard nonfractional differential system 𝛼=𝛼₀=1) to any 𝛼∈(0,𝛼₀] provided that arg(𝜆/𝛼)∈(−𝜋/2,𝜋/2) in the clockwise sense, or equivalently, if |arg𝜆|>𝛼𝜋/2, for all 𝜆∈𝜎(𝐴₀), since ‖‖‖1𝛽₀𝑇⁻¹𝐽_𝐴0𝑇,1𝛽₁𝑇⁻¹𝐴₁𝑇,1𝛽₂𝑇⁻¹𝐴₂𝑇,…,1𝛽_𝑝∀𝛼∈0𝑇,𝛼⁻¹₀𝐴_𝑝𝑇‖‖‖₂<||𝜇₂𝐽_𝐴0𝑑||^1/𝛼0≤||𝜇₂𝐽_𝐴0𝑑||^1/𝛼,.(5.32) Note that the global Lyapunov’s stability conditions (5.23) and (5.26) with nonpositive measures 𝜇₂(𝐽^𝐴1/𝛼_0𝑑) being eventually zero of the corresponding matrices of the unforced fractional dynamic system does not imply the boundedness of the solutions of the system for any admissible forcing bounded control. However, under strict inequalities (5.24) or (5.27) and negative related matrix measures 𝜇₂(𝐽^𝐴1/𝛼_0𝑑), that is, if asymptotic stability holds, the forced solutions for any bounded controls are guaranteed to be uniformly bounded.

It follows after inspecting the solution (3.8), subject to (3.9), and the expressions (5.16) that the stability properties for arbitrary admissible initial conditions or admissible bounded controls are lost in general if 𝛼≥2 and may be improved for 𝛼∈(0,1) compared to the nonfractional calculus counterpart (i.e., for 𝛼=1). However, it turns out that the boundedness of the solutions can be obtained by zeroing some of the functions of initial conditions. Note, in particular, from (5.16) that 𝜑_𝑗is required to be identically zero on its definition domain for 𝑘−1∪{0}∋𝑗<𝛼−1(𝛼≥2) in order that the Γ-functions be positive (note that Γ(𝑥) is discontinuous at zero with an asymptote to −∞ as 𝑥→0⁻). This observation combined with Theorem 5.1 leads to the following direct result which is not a global stability result.

Theorem 5.3. Assume that 𝛼≥2 and the constraint (5.25) holds with negative matrix measure 𝜇₂(𝐽^𝐴1/𝛼_0𝑑). Assume also that 𝜑_𝑗∶[−ℎ,0]→𝐑^𝑛 are any admissible functions of initial conditions for 𝑘−1∪{0}∋𝑗≥𝛼−1 while they are identically zero if  𝑘−1∪{0}∋𝑗<𝛼−1. Then, the unforced solutions are uniformly bounded for all time independent of the delays. Also, the total solutions for admissible bounded controls are also bounded for all time independent of the delays.

The stability of positive or nonnegative solutions is of a direct characterization by combining the positivity conditions of the above section with the stability analysis of this section. The extensions of the given results to discrete fractional systems under either periodic or nonperiodic sampling might be of interest for a future research, [35].

Acknowledgments

The author is grateful to the Spanish Ministry of Education for its partial support of this work through Grant DPI2009-07197. He is also grateful to the Basque Government for its support through Grants GIC07143-IT-269-07, SAIOTEK S-PE08UN15, and SAIOTEK S-PE09UN12.