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

Authors: M. De la Sen
Journal: Abstract and Applied Analysis
 (2011)

---

## 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.


---

## Body

## 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

Z, R,and C are the sets of integer, real, and complex numbers, Z+ and R+ are the positive integer and real numbers,  and
(1.1)Z0+∶=Z+∪{0},R0+∶=R+∪{0},C+∶={z∈C:Rez>0},C0+∶={z∈C:Rez≥0},n̅∶={1,2,…,n}.The  following notation is used to characterize different levels of  positivity of matrices:R0+n×m∶={M=(Mij)∈Rn×m:Mij≥0;∀(i,j)∈n̅×m̅} is the set of all n×m real matrices of nonnegative entries.  If M∈Rn×m then M≥0 is used as a simpler notation for M∈R0+n×m.R+n×m∶={0≠M=(Mij)∈Rn×m:Mij≥0;∀(i,j)∈n̅×m̅} is the set of all nonzero n×m real matrices of nonnegative entries (i.e., at least one of their entries is positive). If M∈Rn×m then M>0 is used as a simpler notation for M∈R+n×m.R++n×m∶={M=(Mij)∈Rn×m:Mij>0;∀(i,j)∈n̅×m̅} is the set of all n×m real matrices of positive entries. If M∈Rn×m then M≫0 is used as a simpler notation for M∈R++n×m. The superscript T denotes the transpose, MiT and Mj are, respectively,  the ith row and the jth column of the matrix M.A close notation to characterize the positivity of vectors is the following:R0+n∶={v=(v1,v2,…,vn)T∈Rn:vi≥0;∀i∈n̅} is the set of all n real vectors of nonnegative components.  If v∈Rn then v≥0 is used as a simpler notation for v∈R0+n. R+n∶={0≠v=(v1,v2,…,vn)T∈Rn:vi≥0;∀i∈n̅} is the set of all  n real nonzero vectors of nonnegative components (i.e., at least one component is positive).  If v∈Rn then v>0 is used as a simpler notation for v∈R+n.R++n∶={v=(v1,v2,…,vn)T∈Rn:vi>0;∀i∈n̅} is the set of all n real vectors of positive components.  If v∈Rn then v≫0 is used as a simpler notation for v∈R++n.M=(Mij)∈Rn×n is a Metzler matrix if Mij≥0;   for all (i,j≠i)∈n̅×n̅.  MRn×n is  the set of Metzler matrices of order n.The maximum real eigenvalue, if any, of a real matrixM, is denoted by λmax(M). 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.


## 1.1. Notation

Z, R,and C are the sets of integer, real, and complex numbers, Z+ and R+ are the positive integer and real numbers,  and
(1.1)Z0+∶=Z+∪{0},R0+∶=R+∪{0},C+∶={z∈C:Rez>0},C0+∶={z∈C:Rez≥0},n̅∶={1,2,…,n}.The  following notation is used to characterize different levels of  positivity of matrices:R0+n×m∶={M=(Mij)∈Rn×m:Mij≥0;∀(i,j)∈n̅×m̅} is the set of all n×m real matrices of nonnegative entries.  If M∈Rn×m then M≥0 is used as a simpler notation for M∈R0+n×m.R+n×m∶={0≠M=(Mij)∈Rn×m:Mij≥0;∀(i,j)∈n̅×m̅} is the set of all nonzero n×m real matrices of nonnegative entries (i.e., at least one of their entries is positive). If M∈Rn×m then M>0 is used as a simpler notation for M∈R+n×m.R++n×m∶={M=(Mij)∈Rn×m:Mij>0;∀(i,j)∈n̅×m̅} is the set of all n×m real matrices of positive entries. If M∈Rn×m then M≫0 is used as a simpler notation for M∈R++n×m. The superscript T denotes the transpose, MiT and Mj are, respectively,  the ith row and the jth column of the matrix M.A close notation to characterize the positivity of vectors is the following:R0+n∶={v=(v1,v2,…,vn)T∈Rn:vi≥0;∀i∈n̅} is the set of all n real vectors of nonnegative components.  If v∈Rn then v≥0 is used as a simpler notation for v∈R0+n. R+n∶={0≠v=(v1,v2,…,vn)T∈Rn:vi≥0;∀i∈n̅} is the set of all  n real nonzero vectors of nonnegative components (i.e., at least one component is positive).  If v∈Rn then v>0 is used as a simpler notation for v∈R+n.R++n∶={v=(v1,v2,…,vn)T∈Rn:vi>0;∀i∈n̅} is the set of all n real vectors of positive components.  If v∈Rn then v≫0 is used as a simpler notation for v∈R++n.M=(Mij)∈Rn×n is a Metzler matrix if Mij≥0;   for all (i,j≠i)∈n̅×n̅.  MRn×n is  the set of Metzler matrices of order n.The maximum real eigenvalue, if any, of a real matrixM, is denoted by λmax(M). 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 thatf:[a,b]→Cn for some real interval [a,b]⊂R satisfiesf∈Ck-2((a,b),Rn) and, furthermore, dk-1f(τ)/dτk-1 exists everywhere in [a,b] for k=[Reα]+1 for some α∈C0+. Then, the Riemann-Liouville left-sided fractional derivative Da+αfRL of order α∈C0+ of the vector function f in [a,b] is pointwise defined in terms of the Riemann-Liouville integral as(2.1)(Da+αfRL)(t)∶=1Γ(k-α)(dkdtk∫atf(τ)(t-τ)α+1-kdτ),t∈[a,b],
where the integer k is given by k=[Reα]+1 and Γ:C∖Z0-→C, where Z0-∶={n∈Z:n≤0}, is the  Γ-function defined by Γ(z)∶=∫0∞τz-1e-τdτ; z∈C∖Z0-. If f∈Ck-1((a,b),Rn) and, furthermore,  f(k)(τ)≡dkf(τ)/dτk exists everywhere in [a,b], then the Caputo left-sided fractional derivative Da+αfC of order α∈C0+ of the vector function f in [a,b] is pointwise defined in terms of the Riemann-Liouville integral as(2.2)(Da+αfC)(t)∶=1Γ(k-α)∫atf(k)(τ)(t-τ)α+1-kdτ,t∈[a,b],
where k=[Reα]+1 if α∉Z0+ and k=α if α∈Z0+. The following relationship between both fractional derivatives holds provided that they exist (i.e., if f:[a,b]→Cn possesses Caputo left-sided fractional derivative in [a,b]), [1](2.3)(Da+αfC)(t)=Da+αRL[f(τ)-∑j=0k-1f(j)(a)(τ-a)jj!](t),t∈[a,b].
Since Reα≤k, the above formula relating both fractional derivatives proves the existence of the  Caputo left-sided fractional derivative in [a,b] if  the Riemann-Liouville one exists in [a,b].


## 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α:(3.1)(D0+αxC)(t)=1Γ(k-α)∫atf(k)(τ)(t-τ)α+1-kdτ=∑i=0pAix(t-hi)+Bu(t)
with k-1<α(∈R+)≤k; k-1, k∈Z0+, 0=h0<h1<h2<⋯<hp=h<∞ being distinct constant delays, A0,Ai∈Rn×n(i∈p̅∶={1,2,…,p}), are the matrices of dynamics for each delay hi,i∈p̅∪{0}, B∈Rn×m is the control matrix. The initial condition is given by k  n-real vector functions φj:[-h,0]→Rn, with j∈k-1¯∪{0}, which are absolutely continuous except eventually in a set of zero measure of [-h,0]⊂R of bounded discontinuities with φj(0)=xj(0)=x(j)(0)=xj0. The function vector u:R0+→Rm is any given bounded piecewise continuous control function. The following result is concerned with the unique solution on R0+ 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 f(t)≡∑i=1pAix(t-hi)+Bu(t)).Theorem 3.1.
The linear and time-invariant differential functional fractional differential system (3.1)  of any order α∈C0+ has a unique solution on R0+ for each given set of initial functions φj:[-h,0]→Rn, j∈k-1¯∪{0} being absolutely continuous except eventually in a set of zero measures of [-h,0]⊂R of bounded discontinuities with φj(0)=xj(0)=x(j)(0)=xj0; j∈k-1¯∪{0} and each given  control  u:R0+→Rm being  a bounded piecewise continuous control function. Such a solution is given by
(3.2)xα(t)=∑j=0k-1(Φαj0(t)xj0+∑i=1p∫0hiΦα(t-τ)Aiφj(τ-hi)dτ)+∑i=1p∫hitΦα(t-τ)Aixα(τ-hi)dτ+∫0tΦα(t-τ)Bu(τ)dτ,t∈R0+
with k=[Reα]+1 if α∉Z+ and k=α if α∈Z+,
(3.3)Φαj0(t)∶=tjEα,j+1(A0tα),Φα(t)∶=tα-1Eα,α(A0tα),(3.4)Eα,j(A0tα)∶=∑l=0∞(A0tα)lΓ(αl+j),j∈k-1¯∪{0,α}
for t≥0 and Φα0(t)=Φα(t)=0 for t<0, where Eα,j(A0tα) 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: (3.5)(D0+αxRL)(t)=∑i=0pAix(t-hi)+Bu(t)
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(3.6)(D0+αxC)(t)=∑i=0pAix(t-hi)+Bu(t)-D0+αRL(∑j=0k-1xj0τjj!)(t).
Since the Caputo left-sided fractional derivative and the Riemann-Liouville fractional integral of order α∈Ĉ+∶={Z+∪{z∈C+:Rez∉Z+}} are inverse operators (what is not the case if α∉Ĉ+), (see [1, Lemma 2.21(a)]), one gets from (3.6), (2.3),  and (3.2)  if α∈Ĉ the subsequent  result for the fractional differential system (3.5) on R0+.Corollary 3.2.
If (3.5) of any order α∈Ĉ+ is replaced with (3.1) under the same initial conditions then its unique solution on R0+ is given by
(3.7)xα(t)=∑j=0k-1((Φαj0(t)-tjj!In)xj0+∑i=1p∫0hiΦα(t-τ)Aiφj(τ-hi)dτ)+∑i=1p∫hitΦα(t-τ)Aixα(τ-hi)dτ+∫0tΦα(t-τ)Bu(τ)dτ,α∉Z0+;t∈R0+=∑j=0k-1((Eα,j+1(t)-1j!In)tjxj0+∑i=1p∫0hiΦα(t-τ)Aiφj(τ-hi)dτ)+∑i=1p∫hitΦα(t-τ)Aixα(τ-hi)dτ+∫0tΦα(t-τ)Bu(τ)dτ,α∉Z0+,t∈R0+
with k=[Reα]+1 if α∉Z+ and k=α if α∈Z+.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 (f(t)≡Bu(t)) 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:
(3.8)xα(t)=∑j=0k-1(Ψαj0(t)xj0+∑i=1p∫0hiΨαj0(t-τ)φj(τ-hi)dτ)+∫0tΨα(t-τ)Bu(τ)dτ
for  t∈R0+, any α∈C+ with k=[Reα]+1 if α∉Z+ and k=α if α∈Z+;
(3.9)Ψαj0(t)∶=tjEα,j+1(A0tα)+∑i=1p∫0tτα-1Eαα(A0τα)AiΨαj0(t-τ-hi)dτ,Ψα(t)∶=tα-1Eαα(A0tα)+∑i=1p∫0tτα-1Eαα(A0τα)AiΨα(t-τ-hi)dτ
for t≥0 and Ψαj0(t)=Ψα(t)=0, j∈k-1¯∪{0} for t∈[-h,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 k=[Reα]+1 if α∉Z+ based on (3.6).Corollary 3.4.
If (3.5) being of order α∈Ĉ+ is replaced with (3.1) under the same initial conditions then its unique solution on R0+ is given by(3.10)xα(t)=∑j=0k-1((Ψαj0(t)-tjj!In)xj0+∑i=1p∫0hiΨαj0(t-τ)φj(τ-hi)dτ)+∫0tΨα(t-τ)Bu(τ)dτ
with k=[Reα]+1 if α∉Z+ and k=α if α∈Z+ 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)
α=k which yields the solution:(3.11)xk(t)=∑j=0k-1(Ψkj0(t)xj0+∑i=1p∫0hiΨkj0(t-τ)φj(τ-hi)dτ)+∫0tΨk(t-τ)Bu(τ)dτ,(2)
a further particular caseα=k=1 which yields the solution:(3.12)x1(t)=Ψ1(t)xj0+∑i=1p∫0hiΨ1(t-τ)φ0(τ-hi)dτ+∫0tΨ1(t-τ)Bu(τ)dτ
since Ψ100(t)=Ψ1(t), t∈R0+ which is the unique solution of (Dx)(t)=∑i=0pAi(t-hi)+Bu(t) under any almost everywhere absolutely continuous function (except eventually in some subset of zero measure of [-h,0] of bounded discontinuities) of initial conditions  φ≡φ0:[-h,0]→Rn. Use for this case, the less involved notations (Ψx)(t)=(Ψ100x)(t)=(Ψ1x)(t) for the smooth evolution operator from R0+×Rn to Rn, and Φ(t)=Φ100(t)=Φ1(t)=eA0t, t∈R+  for the exponential matrix function eA0t from R0+ to Rn×n, which defines a C0-semigroup (eA0t,t∈R0+) of infinitesimal generator A0 from R0+ to L(Rn). Then, the unique solution x(t)≡x1(t),t∈R+ for the given function of initial conditions is(3.13)x(t)=Φ(t)x0+∑i=1p∫0hiΦ(t-τ)Aiφ(τ-hi)+∑i=1p∫hitΦ(t-τ)Aix(τ-hi)dτ+∫0tΦ(t-τ)Bu(τ)dτ=Ψ(t)x0+∑i=1p∫0hiΨ(t-τ)φ(τ-hi)dτ+∫0tΨ(t-τ)Bu(τ)dτ,t∈R0+
and x(t)=φ(t) for t∈[-h,0], where Φ(t)=eA0t satisfies Φ̇(t)=A0Φ(t)t∈R, and Ψ̇(t)=∑i=0pAiΨ(t-hi) with Ψ(0)=Φ(0)=In (the n-identity matrix) and Ψ(t)=0, t∈[-h,0) which has a unique solution Ψ(t)=eA0t(In+∑i=1p∫hite-A0τAiΨ(τ-hi)dτ)  for t∈R0+, [2, 3, 25, 26].A problem of interest  when considering a set ofp delays in [0,h] is the case of  potentially repeated delays, then subject to 0=h0≤h1≤h2≤⋯≤hp≤h<∞, with q of them {hjp,j∈q¯∪{0}} being distinct, each being repeated 1≤νj≤p(j∈q¯∪{0}) times  so that(3.14)0=h0=h0p<h1p<⋯<hqp≤h<∞,∑j=0pνjp=p+1,hjp=hk+i,k=∑l=0j-1νl,∀i∈ν¯j,j∈q¯∪{0}.Thus, the following result holds from Theorem3.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 p≥0 (potentially repeated) delays and 0≤q≤p distinct  delays:(3.15)(D0+αxC)(t)=∑i=0pAi(t-hi)+Bu(t),(D0+αxRL)(t)=∑i=0pAi(t-hi)+Bu(t)
on R0+ for the given set of initial conditions on [-h,0) are given by
(3.16)xα(t)=∑j=0k-1(Ψαj0(t)xj0+∑i=1q∫0hipΨαj0(t-τ)φj(τ-hip)dτ)+∫0tΨα(t-τ)Bu(τ)dτ
for any α∈C+ with k=[Reα]+1 if α∉Z0+ and k=α if α∈Z0+, and, respectively by
(3.17)xα(t)=∑j=0k-1((Ψαj0(t)-tjj!)xj0+∑i=1q∫0hipΨαj0(t-τ)φj(τ-hip)dτ)+∫0tΨα(t-τ)Bu(τ)dτ
for any α∈Ĉ+ with k=[Reα]+1 if α∉Z+ and k=α if α∈Z+, where
(3.18)Ψαj0(t):=Eα,j+1((∑i=0ν0-1Ai)tα)+∑i=1q∫0tτα-1Eα,α((∑i=0ν0-1Ai)τα)(∑l=1νiA(∑j=0i-1νj+l))Ψαj0(t-τ-hi)dτ,t≥0 and Ψαj0(t)=Ψα(t)=0, j∈k-1¯∪{0} for t∈[-h,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 φj:[-h,0]→R0+n, j∈k-1¯∪{0}.Theorem 4.1.
The Caputo fractional differential  system (3.1) under the delay constraint 0=h0<h1<h2<⋯<hp=h<∞ for any  given absolutely continuous functions of initial conditions φj:[-h,0]→R0+n, j∈k-1¯∪{0} and any piecewise continuous vector function u:R0+→R0+n if k=[α]+1 if α∉Z+ and k=α∈Z+,forallt∈R0+; for all α∈R+ has  following properties:
(i)
Φαj0(t) is nonsingular; forallj∈k-1¯∪{0} and Φα(t)≥0; for all  t∈R0+ (if B∈R+n×m then Φα(t)>0; for all t∈R0+),(ii)
(1)
A0∈MRn×n⇔Φ(t)≡Φ100(t)>0; for all t∈R0+,(2)
A0∈MRn×n⇒Φαj0(t)>0;forallj∈k-1¯∪{0}; for all t∈[0,t̅j) for some sufficiently small t̅j∈R+ with Φα00(t)>0, for all  t∈R0+(i.e.,t̅0=∞). This property holds for all  t∈R0+  (i.e., t̅j=∞; for all j∈k-1¯∪{0}) if, in addition, either A0≥0 or if A0 is nilpotent or if 0<α≤k=1. Furthermore, there are at least n entries (one per row) of Φαj0(t) being positive; forallt∈R0+;(iii)
Any solution (3.2) to any Caputo fractional differential system (3.1) is nonnegative independent of the delays; that is, xα(t)∈R0+n; for all t∈[-h,t̅)∩R0+ for some t̅∈R0+, for any set of delays satisfying  0=h0<h1<h2<⋯<hp≤h<∞ and any absolutely continuous functions of initial conditions φj:[-h,0]→R0+n,  for all  j∈k-1¯∪{0} and any piecewise continuous control u:R0+→R0+m, if and only if A0∈MRn×n for t̅∈R0+ being sufficiently small. Furthermore, xα(t)∈R0+n; for all  t∈[-h,0)∪R0+ if, in addition,  either A0≥0 or if A0 is nilpotent or if  0<α≤k=1,Ai∈R0+n×n(∀i∈p̅) and B∈R0+n×m.Proof.
It is now proven thatΦαj0(t)≥0; for all t∈R0+⇒A0∈MRn×n; for all  α∈R+ for any j∈k-1¯∪{0}. First, note the following. If α=k=1 then Φα00(t)=E1,1(A0t)=Φ(t)=∑ℓ=0∞A0ℓtℓ/ℓ!=eA0t≥0 if A0∈MRn×n from the above part of the proof and also A0∈MRn×n⇒Φ(t)≥0; for all t∈R0+. This follows by contradiction. Assume that Φim(t)<0 for some t∈R+. Consider the positive differential system ẋ(t)=A0x(t), x(0)=ej, A0∈MRn×n so that xi(t)=-|Φim(t)|<0 which contradicts the system being positive. Thus, A0∈MRn×n⇔Φ(t)≥0; for all  t∈R0+. Furthermore,  since Φ(t) 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:
A0∈MRn×n⇔(Φ(t)=Φα00(t)=eA0t>0∧Φ(t) is non-singular; for all t∈R0+). Since Φ(t) is nonsingular; for all  t∈R0+ at least n of its entries (one per-row) is positive. Property (i) has been proven. Now, one gets from (3.3)-(3.4):
(4.1)Φαj0(t)=Eα,j+1(A0tα)=∑l=0∞A0ltαlΓ(αl+j+1)=∑l=0∞A0ltll!t(α-1)ll!(αl+j)Γ(αl+j),j∈k-1¯∪{0}.
Let ei the ith unit Euclidean vector of  Rn whose ith component is 1. Then, one obtains for all j∈k-1¯∪{0}, irrespective of the value of  α∈R+ and k∈Z+ being k=[α]+1 if α∉Z+ and k=α∈Z+, provided that A0∈MRn×n:
(4.2)(Φαj0im(t))=eiTΦαj0(t)em=eiT(∑l=0∞A0ltll!t(α-1)ll!Γ(αl+j+1))em=eiT(∑l=0∞A0ltll!t(α-1)ll!(αl+j)Γ(αl+j))em≥eiT(eA0t)emmin⁡0≤l≤∞(t(α-1)ll!Γ(αl+k))=eiT(∑l=0∞A0ltll!t(α-1)ll!(αl+j)Γ(αl+j))em(4.3)=eiT(∑l=0NA0ltll!)emmin⁡0≤l≤N(t(α-1)ll!(αl+j)Γ(αl+k))=eiT(eA0t)emmin⁡0≤l≤N(t(α-1)lΓ(l+1)(αl+j)Γ(αl+k))≥eiT(eA0t)emmin⁡0≤l≤N(t(α-1)lΓ(l+1)(αl+j)Γ(αl+k))
for all (i,m)∈n̅×n̅, for all t∈R0+ since A0∈MRn×n⇔Φ(t)≥0; for all t∈R0+, for some N(≤∞)∈Z0+ and N is finite if and only if A0 is nilpotent (of degree N). Equation (4.3) implies that eiT(eA0t)em≥0 and then (Φαj0im(t))≥0, for all (i,m)∈n̅×n̅ in the following cases:
(a)N≤∞, t∈[0,t̅), since ∩R0+ and some sufficiently small t̅∈R+,  since A0∈MRn×n⇒I+A0t>0; for all  t∈R0+ and Φαj0(t)=∑ℓ=0N(A0ℓtℓ/ℓ!)(t(α-1)ℓℓ!/(αℓ+j)Γ(αℓ+j))=I+A0t+o(t)>0 for some sufficiently small t∈R0+ for any α∈R+.
(b)N≤∞ and A0≥0 since  t(α-1)ℓℓ!/(αℓ+j)Γ(αℓ+j)≥0; forallℓ∈Z0+, for all  j∈k-1¯∪{0} for any α∈R+. It follows from inspection of (4.2) since eiT(eA0t)em≥0 for all (i,m)∈n̅×n̅, since A0∈MRn×n. This implies (Φαj0im(t))≥0;  for all t∈R0+.
(c)N<∞⇒min0≤ℓ≤N(t(α-1)ℓ/Γ(αℓ+k))>0; for all  j∈k-1¯∪{0} for any α∈R+  so that eiT(eA0t)em≥0, for all (i,m)∈n̅×n̅, since A0∈MRn×n,  irrespectively of A0≥0 or not, what follows from (4.3). This implies (Φαj0im(t))≥0; for all  t∈R0+.
(d)N≤∞, 0<α≤k=1. Then, j=0 so that
(4.4)t(α-1)ll!Γ(αl+j+1)=t(α-1)ll!Γ(αl+1)=t(α-1)ll!αlΓ(αl)=t(α-1)l(l-1)!αΓ(αl)=Γ(l)αt(1-α)lΓ(αl)≥1αt(1-α)l,∀l∈Z0+,∀t∈R0+
since 0<α≤1 implies
(4.5)Γ(αl)=∫0∞ταl-1e-τdτ≤Γ(l)=∫0∞τl-1e-τdτ,∀l∈Z0+.
As a result, Φα00(t) from (4.2); for all  t∈R0+. Also, direct calculations with (3.3)-(3.4) lead to
(4.6)Φα(t)∶=tα-1Eαα(A0tα)=∑l=0∞tα-1A0ltαlΓ((l+1)α)=∑l=0∞A0ltll!t(α-1)(l+1)l!Γ((l+1)α)
and similar developments to the above ones yield (Φα(t))im≥0; for all  (i,m)∈n̅×n̅, for all t∈R0+ under the same conditions as above in the cases (a) to (d) for Φα00(t). On the other hand, one gets from (3.2)–(3.4) for the unforced system with point initial conditions at t=0:
(4.7)xα(t)=∑j=0k-1Φαj0(t)xj0=[Φα00(t),…,Φα,k-1,0(t)][x00T,…,xk-1,0T]T
which leads to xα(t)=Φαi0(t)xi0 by taking point initial conditions xi0≠0,xj0=0,(i≠j),   j∈k-1¯∪{0} so that Φαi0(t) is nonsingular for all t∈R0+ 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 xi0≠0, xj0=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 j∈k-1¯∪{0}, Φαj0(t) is nonsingular for all  j∈k-1¯∪{0}, all  and, in addition, Φαj0(t)>0; forallj∈k-1¯∪{0}, for all t∈R0+ if either A0≥0 or if A0 is nilpotent or if  0<α≤k=1 or without these restricting condition within  some first interval [0,t̅). The following properties have been proven:
(a)A0∈MRn×n⇔Φα00(t)>0; for all  t∈R0+,
(b)A0∈MRn×n⇒(Φαj0(t)>0∧detΦαj0(t)≠0∧Φα(t)≥0; for all  j∈k-1¯∪{0}, k=[α]+1 if α∉Z+and k=α∈Z+,  for all t∈R0+).
It remains to proveΦαj0(t)>0; for all t∈R0+⇒A0∈MRn×n; for all  j∈k-1¯, some t∈R0+. This is equivalent to its contrapositive logic proposition. Proceed by contradiction by assuming thereexistj∈k-1¯ such that A0∉MRn×n⇒Φαj0(t)<0, some t∈R0+. Note that  A0∉MRn×n⇒eiTeA0t(t)em=eiTΦα00(t)em<0, some t∈R0+,  some (i,m≠i)∈n̅.  Then, one gets
(4.8)(Φαj0im(t))=eiTΦαj0(t)em≤eiT(eA0t)emmax⁡0≤l≤∞(t(α-1)ll!Γ(αl+k))<0
which contradicts (-Φαj0(t))>0; for allt∈R0+⇐A0∉MRn×n; for all  j∈k-1¯, some t∈R0+. Thus, the proof of Properties (i)-(ii) becomes complete since the above proven  property (a) extends to any j∈k-1¯∪{0} as follows.
(c)A0∈MRn×n⇔Φαj0(t)>0; for all  j∈k-1¯∪{0}, k=[α]+1 if α∉Z+ and k=α∈Z+, forallt∈R0+; for all α∈R+so that the unforced solution for any set of nonnegative point initial conditions is nonnegative for all time and, furthermore, φi(t)≥0(∀i∈k-1¯∪{0}); for all t∈[-h,0],u(t)∈R0+n; for all t∈R0+,Ai≥0(∀i∈p̅)andB≥0; for all t∈R0+ 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 Ai(somei∈p̅)  which is negative, or if A0∉MRn×n, it can always be found a control u(t)∈R0+n 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 t=0 such that some component of the solution is negative at t=0+. Thus, Property (iii) is proven.The following result is obvious from the proof of Theorem4.1.Corollary 4.2.
Theorem4.1(iii) is satisfied also independent of the delays for any given  set of delays satisfying the constraint 0=h0≤h1≤h2≤⋯≤hp=h<∞.Proof.
It follows directly since Theorem4.1 is an independent of the delay size type result  and, under the delay constraint 0=h0≤h1≤h2≤⋯≤hp=h<∞, it has also to be fulfilled for any combination of delays satisfying the  stronger constraint   0=h0<h1<h2<⋯<hp=h<∞.Corollary 4.3.
Any solution (3.8), subject to (3.9),  to the Caputo fractional differential  system (3.1) under the delay constraint 0=h0≤h1≤h2≤⋯≤hp=h<∞ is nonnegatively independent of the delays  within a first interval, that is, it satisfies  xα(t)∈R0+n; for all   t∈[-h,t̅)∩R0+ for some sufficiently small t̅∈R0+ for any given absolutely continuous functions of initial conditions φj:[-h,0]→R0+n, j∈k-1¯∪{0}  and any given piecewise continuous vector function u:R0+→R0+n with k=[α]+1 if α∉Z+ and k=α∈Z+, for all  t∈R0+; for all  α∈R+ if and only if A0∈MRn×n, Ai∈R0+n×n(∀i∈p̅), and B∈R0+n×m. In addition,  xα:[-h,0)∪R0+→R0+n if, in addition,  either A0≥0 or if A0 is nilpotent or if 0<α≤k=1. Furthermore, Ψαj0(t)>0 (with at least n entries being positive), detΨαj0(t)>0(∀j∈k-1¯∪{0}) and Ψα(t)≥0; for allt∈R0+(ifB∈R+n×mthenΨα(t)>0;∀t∈R0+).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 somet∈R+ provided that B≫0 and the control u:R0+→R0+n is admissible (i.e., piecewise continuous) and nonidentically zero since Ψα(t)>0 and nonsingular for all t∈R+. 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 φj:[-h,0]→R0+n, j∈k-1¯∪{0} and any given piecewise continuous vector function u:R0+→R0+n with k=[α]+1 if α∉Z+ and k=α∈Z+, for all  t∈R0+; for all α∈R+. 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=h0≤h1≤h2≤⋯≤hp=h<∞ is everywhere nonnegative independent of the delays, that is, it satisfies xα:[-h,0)∪R0+→R0+n, for any given absolutely continuous functions of initial conditions φj:[-h,0]→R0+n, j∈k-1¯∪{0} and any given piecewise continuous vector function u:R0+→R0+n with k=[α]+1 if α∉Z+ and k=α∈Z+, forallt∈R0+; for all α∈R+ if A0∈MRn×n, Ai∈R0+n×n(∀i∈p̅), (Eα,j+1(t)-(1/j!)In)≥0; forallj∈k-1¯∪{0}, for all  t∈R0+ and B∈R0+n×m. The conditions A0∈MRn×n, Ai∈R0+n×n(∀i∈p̅) and B∈R0+n×m are also necessary for  xα:[-h,0)∪R0+→R0+n for any nonnegative function of initial conditions and nonnegative controls. The condition (Eα,j+1(t)-(1/j!)In)≥0; for all j∈k-1¯∪{0}, for all t∈R0+ is removed for initial conditions φj:[-h,0]→R0+n subject to φj(0)=xj0=0.Proof.
The proof follows  in a similar way as the sufficiency part of the proof of  Theorem4.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=h0≤h1≤h2≤⋯≤hp=h<∞ is everywhere nonnegatively independent of the delays, that is, it satisfies xα:[-h,0)∪R0+→R0+n, for any given absolutely continuous functions of initial conditions φj:[-h,0]→R0+n, j∈k-1¯∪{0} and any given piecewise continuous vector function u:R0+→R0+n with k=[α]+1 if α∉Z+ and k=α∈Z+, for all t∈R0+; forallα∈R+ if and only if A0∈MRn×n, Ai∈R0+n×n(∀i∈p̅), (Ψαj0(t)-(tj/j!)In)≥0; for all j∈k-1¯∪{0}, for all t∈R0+ and B∈R0+n×m. The condition (Ψαj0(t)-(tj/j!)In)≥0; for all j∈k-1¯∪{0}, for all t∈R0+ is removed for initial conditions φj:[-h,0]→R0+n subject to φj(0)=xj0=0.Proof.
The proof of sufficiency follows  in a similar way as the sufficiency part of the proof of  Theorem4.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 thatA0∉MRn×n and the solution is nonnegative for all time for any nonnegative function of initial conditions and controls. Take initial conditions φj(t)=0; for all t∈[-h,0], for all j∈k-1¯; φ0(t)=0; for all t∈[-h,0), φ0(0)=x00(0)≠0and u≡0 on R0+. Then (3.10) becomes
(4.9)xα(t)=(Ψα00(t)-In)x00=(Eα00(t)-In)x00fort∈[0,h1]
since Ψα00(t)=Φα00(t) for t∈[0,h1]. Since A0∉MRn×n, there exist  t∈[0,h1] and (i=i(t),m(t)≠i)∈n̅ such that (Φα00(t))im<0. Otherwise, if A0∉MRn×n and Φα00(t)>0; for all t∈[0,h1], it would follow from (4.3) that Φα00(t)>0; for all t∈R0+ since
(4.10)eA0t=eχA0h1+A0δ=eχA0h1eA0δ=(eA0h1)χeA0δ>0
from the semigroup property of (eA0t,t∈R0+) with χ=χ(t,h1)=[t/h1] and (0,h1)∋δ=δ(t,h1)=t-χh1 what implies Φα00(t)>0; for all t∈R0+ from (4.3). Thus, A0∈MRn×n which contradicts A0∉MRn×n. It has been proven that  A0∉MRn×n⇒eiTΦα00(t)em<0; forallt∈(0,h1] for some (i=i(t),m(t)≠i)∈n̅. Now, take x00j=δjm(∀j∈n̅) where δjm denotes the Kronecker delta. Then,
(4.11)xαi(t)=eiT(Ψα00(t)-tjj!In)x00=eiT(Φα00(t)-tjj!In)em=eiTΦα00(t)em<0.
As a result, A0∈MRn×n 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 thateiTAℓej<0 and  hℓ≠hi; for all i(≠ℓ)∈p̅ for some i,j∈n̅, ℓ∈p̅. Take initial conditions  xj0=φj(0)=0; forallt∈[-h,0]; forallj∈k-1¯∪{0}, φj≡0; forallj∈k-1¯ and  u≡0. One gets from (3.2)
(4.12)xαi(t)=∫0h′eiTΦα(t-τ)(∑j(≠l)∈p̅Aj)φ0(τ-h′)dτ+∫0hleiTΦα(t-τ)Alφ0(τ-hl)dτ;∀t∈[0,hi]
for the case h′=hi;  foralli(≠ℓ∈p̅). Now, if h=hℓ>h′, take a further specification of initial conditions as follows: φ0(t)=0; forallt∈[0,h′], and φ0(τ)=(k1,…,kn)T≫0; forallt∈(h′,hℓ] then
(4.13)xαi(t)=∫h′hleiTΦα(t-τ)Alφ0(τ-hl)dτ=∑r=1n∑m=1n(∫h′hlΦαir(t-τ)Alrmdτ)km=∑j=1nkj(∫h′hlΦαiT(t-τ)dτ)Alj=∑(m≠j)∈n̅km(∑r=1m∫h′hlΦαir(t-τ)dτ)Alrm+kj(∑r=1m∫h′hlΦαir(t-τ)dτ)Alrj,∀t∈[0,hi].
As a result,   Ai≥0(∀i∈p̅) 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, andB≥0 is not fulfilled so that it exists at least an  entry Bℓj<0 of B. Then, one has under identically zero initial conditions the following unique solution:
(4.14)xαi(t)=∫0teiTΨα(t-τ)Bu(τ)dτ=∑i=1m∫0tΨαiT(t-τ)Biui(τ)dτ=∑i=1m∑l=1n∫0tΨαil(t-τ)Bliui(τ)dτ=∑(i≠j)∈m̅∑l=1n∫0tΨαil(t-τ)Bliui(τ)dτ-(∑l=1n∫0tΨαjl(t-τ)|Blj|dτ)kuj<0,t∈R0+
provided that kuj>∑(i≠j)∈m̅∑ℓ=1n∫0tΨαiℓ(t-τ)Bℓiui(τ)dτ/(∑ℓ=1n∫0tΨαjℓ(t-τ)|Bℓj|dτ) by assuming that B≥0 fails because Bℓj<0 for some (ℓ,j)∈n̅×m̅ and a constant control component uj≡kuj>0 is injected on the time interval [0,t] for some arbitrary t∈R+ for the remaining control components being chosen t  be nonnegative for all time. This contradicts that the solution is nonnegative for all time if the condition B≥0 fails.Remark 4.6.
Note that Theorem4.1 can be extended as a necessary condition for t∈[0,h1]  since Ψαj0(t)=Φαj0(t) for t∈[0,h1]; forallj∈k-1¯∪{0}, forallt∈R0+.Remark 4.7.
Note  by simple calculation that(A0∈MRn×n∧Ai≥0(i∈p̅))⇒(∑i=0pAi)∈MRn×n. This is a necessary and sufficient condition for the nonnegativity of the solutions of the Caputo fractional differential system (3.1) of arbitrary order α∈R+ under arbitrary nonnegative controls and initial conditions  in the absence of delays; that is, for hi=0; foralli∈ω̅∪{0} and any ω∈Z+.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 ofp delays. However, the conditions are weakened for particular situations involving repeated delays as follows. Note from Theorem 4.5 that the various given conditions Ai≥0 of necessary type to guarantee the nonnegativity of the solution under any admissible nonnegative controls and nonnegative initial conditions are weakened to (∑ℓ=1νiA(∑j=0i-1νj+ℓ)) if there is some repeated delay hi of multiplicity νi≥2 (i.e., the number of distinct delays is 0≤q<p). Also, if h0=0 is repeated with multiplicity ν0≥2 then the condition A0∈MRn×n for ν0=1 is replaced by (∑ℓ=0ν0-1Aℓ)∈MRn×n.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α∈R+

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 α∈R+, (see (1.8.27)–(1.8.29), [1]).(1)
If0<α<2 then for |z|→∞ and some μ∈R satisfying μ<πmin(1,α):(5.1)Eαβ(z)=1αz(1-β)/αe(z1/α)-∑j=1N1Γ(β-αj)1zj+O(1zN+1)
with |argz|≤μ<πmin(1,α), any N∈Z+, and (5.2)Eαβ(z)=-∑j=1N1Γ(β-αj)1zj+O(1zN+1)
with π≥|argz|≥μ<π=πmin(1,α), any N∈Z+.(2)
Ifα≥2 then for |z|→∞(5.3)Eαβ(z)=1α∑j∈Q(z1/αe2jπi/α)1-βe(e2jπi/αz1/α)-∑j=1N1Γ(β-αj)1zj+O(1zN+1)
for any N∈Z+ with β∈k-1¯∪{0}∪{α}, |argz|≤απ/2, Q∶={n∈Z:|argz+2πn|≤απ/2} and i=-1 being the complex imaginary unit. The above formulas are extendable to the Mittag-Leffler matrix functions Eα,j+1(A0tα)∶=∑ℓ=0∞(A0tα)ℓ/Γ(αℓ+j+1); forallj∈k-1¯∪{0}, respectively, Eαα(A0tα) by identifying z→A0tα, z1/α→(A0)1/αt (if (A0)1/αexists) and z-1→A0-1t-α (if A0 is non-singular), β→j+1, respectively,  β→α.  Irrespective of the existence of (A0)1/αand of A0≠0 being singular or nonsingular, it is possible to identify z-1→∥A0∥-1t-α and z→∥A0∥tα and to use(5.4)‖Eα,j+1(A0tα)‖≤Eα,j+1(‖A0‖tα)=∑l=0∞(‖A0‖tα)lΓ(αl+j+1),∀j∈k-1¯∪{0},‖Eαα(A0tα)‖≤Eαα(‖A0‖tα)=∑l=0∞(‖A0‖tα)lΓ(αl+α).The method may be used to calculate an asymptotic estimate of the solution (3.2) if A0 is non-singular (or an upperbounding function for any nonzero A0) 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 A0 appear by conjugate pairs A0 then A0=T-1J0A0 where J0 is its real canonical form. First, consider two separate cases as follows.(A) Assume thatα∈R+, A0 is real non-singular and (A0)1/α exists; that is, there exist M such that Mα=A0 and Ai(i∈p̅) is real. Then, one gets from (5.1)–(5.3):(5.5)Eα,j+1(A0tα)=1α(A0-j)1/αt-je(A01/α)t-∑l=1N1Γ((1-α)l+1)(A0-l)t-lα+O(A0-(N+1)t-(N+1)α),∀j∈k-1¯∪{0},Eαα(A0tα)=1αA0(1-α)/αt1-αe(A01/α)t-∑l=1N1Γ((1-l)α)(A0-l)t-lα+O(A0-(N+1)t-(N+1)α)
as t→∞ if 0<α<2, for any N∈Z+, (5.6)Eα,j+1(A0tα)=1α∑l∈Q(A01/αte2lπi/α)-je(e2lπi/α(A01/α)t)-∑l=1N1Γ((1-α)l+1)(A0-l)t-lα+O(A0-(N+1)t-(N+1)α),∀j∈k-1¯∪{0},Eαα(A0tα)=1α∑l∈Q(A01/αte2lπi/α)1-αe(e2lπi/α(A01/α)t)-∑l=1N1Γ((1-l)α)(A0-l)t-lα+O(A0-(N+1)t-(N+1)α)
as t→∞ if α≥2, for any N∈Z+, with Q∶={n∈Z:|n|≤α/4}.(B) Assume thatα∈R+ and Ai(i∈p̅∪{0}) is real, one obtains from (5.1)-(5.2):(5.7)|Eαβ(z)|≤E̅αβ(z)=1α|z|(1-β)/α|ez|1/α+|∑l=1N1Γ(β-αl)1zj+O(1zN+1)|
for 0<α<2, β∈k-1¯∪{0},(5.8)|Eαβ(z)|≤E̅αβ(z)=1α∑j∈Q(|z|(1-β)/α)|ez|1/α+|∑l=1N1Γ(β-αl)1zj+O(1zN+1)|
for α≥2, β∈k-1¯∪{0,α}. Thus, on gets from (5.7)(5.9)‖Eα,j+1(A0tα)‖≤1α(‖A0‖-j)1/αt-j‖eA0t‖1/α-∑l=1N1Γ((1-α)l+1)(‖A0‖-l)t-lα+O(‖A0‖-(N+1)t-(N+1)α),∀j∈k-1¯∪{0},‖Eαα(A0tα)‖≤1α‖A0‖(1-α)/αt1-α‖eA0t‖1/α-∑l=1N1Γ((1-l)α)(‖A0‖-l)t-lα+O(‖A0‖-(N+1)t-(N+1)α)
as t→∞, for any N∈Z+, if 0<α<2, and one gets from (5.8)(5.10)‖Eα,j+1(A0tα)‖≤1α∑l∈Q(‖A0‖1/αt)-j‖eA0t‖1/α-∑l=1N1Γ((1-α)l+1)(‖A0‖-l)t-lα+O(‖A0‖-(N+1)t-(N+1)α),∀j∈k-1¯∪{0},‖Eαα(A0tα)‖≤1α∑l∈Q(‖A0‖1/αt)1-α‖eA0t‖1/α-∑l=1N1Γ((1-l)α)(‖A0‖-l)t-lα+O(‖A0‖-(N+1)t-(N+1)α)
as t→∞ if α≥2, for any N∈Z+,  with Q∶={n∈Z:|n|≤α/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α∈R+, A0 is real non-singular, (A0)1/α exists and Ai(i∈p̅) are also real. Then, one gets from (5.5)–(5.6) into (3.9):(5.11)Ψαj0(t)=1α(A0-j)1/αe(A01/α)t-∑l=1N1Γ(j+1-αl)(A0-l)tj-lα+O(A0-(N+1)tj-(N+1)α)+∑i=1p∫0t(1αA0(1-α)/αe(A01/α)τ-∑l=1N1Γ((1-l)α)(A0-l)τ(1-l)α-1+O(A0-(N+1)τ-Nα-1))×AiΨαj0(t-τ-hi)dτ,(5.12)Ψα(t)=1αA0(1-α)/αe(A01/α)t-∑l=1N1Γ((1-l)α)(A0-l)t(1-l)α+O(A0-(N+1)t-Nα)+∑i=1p∫0t(1αA0(1-α)/αe(A01/α)τ-∑l=1N1Γ((1-l)α)(A0-l)τ(1-l)α-1+O(A0-(N+1)τ-Nα-1))×AiΨα(t-τ-hi)dτ
for all j∈k-1¯∪{0} as t→∞ if 0<α<2, for any N∈Z+, and(5.13)Ψαj0(t)=1α∑l∈Q(A01/αe2lπi/α)-je(e2lπi/α(A01/α)t)-∑l=1N1Γ(j+1-αl)(A0-l)tj-lα+O(A0-(N+1)tj-(N+1)α)+∑i=1p∫0t(1α∑l∈Q(A0(1-α)αe2lπiα)-je(e2lπiα(A01α)τ)1α∑l∈Q(A0(1-α)/αe2lπi/α)-je(e2lπi/α(A01/α)τ)-∑l=1N1Γ((1-l)α)(A0-l)τ(1-l)α-1+O(A0-(N+1)τ-Nα-1)1α∑l∈Q(A0(1-α)αe2lπiα)-je(e2lπiα(A01α)τ))×AiΨαj0(t-τ-hi)dτ,Ψα(t)=1α∑l∈Q(A0(1-α)/αe2lπ(1-α)i/α)e(e2lπi/α(A01/α)t)-∑l=1N1Γ((1-α)l+1)(A0-l)tj-lα+O(A0-(N+1)tj-(N+1)α)+∑i=1p∫0t(1α∑l∈Q(A0(1-α)αe2lπ(1-α)iα)e(e2lπiα(A01α)τ)1α∑l∈Q(A0(1-α)/αe2lπ(1-α)i/α)e(e2lπi/α(A01/α)τ)-∑l=1N1Γ((1-α)l)(A0-l)τ(1-l)α-1+O(A0-(N+1)τ-Nα-1)1α∑l∈Q(A0(1-α)αe2lπ(1-α)iα)e(e2lπiα(A01α)τ))×AiΨαj0(t-τ-hi)dτ
for all j∈k-1¯∪{0} as t→∞ if α≥2, for any N∈Z+.(2)  Assume thatα∈R+ and Ai(i∈p̅∪{0}) are real. Then, (5.14)‖Ψαj0(t)‖≤1α(‖A0‖-j)1/α‖eA0t‖1/α+∑l=1N1|Γ(j+1-αl)|‖A0‖-ltj-lα+O(‖A0‖-(N+1)tj-(N+1)α)+∑i=1p∫0t(1α‖A0‖(1-α)α‖eA0τ‖1α+∑l=1N1|Γ((1-l)α)|(‖A0‖-l)τ(1-l)α-11α‖A0‖(1-α)/α‖eA0τ‖1/α+∑l=1N1|Γ((1-l)α)|(‖A0‖-l)τ(1-l)α-1+O(‖A0‖-(N+1)τ-Nα-1)1α‖A0‖(1-α)α‖eA0τ‖1α+∑l=1N1|Γ((1-l)α)|(‖A0‖-l)τ(1-l)α-1)‖Ai‖‖Ψαj0(t-τ-hi)‖dτ,(5.15)‖Ψα(t)‖≤1α‖A0‖(1-α)/α‖eA0t‖1/αtα+∑l=1N1|Γ((1-l)α)|(‖A0‖-l)t(1-l)α+O(‖A0‖-(N+1)t-Nα)+∑i=1p∫0t(1α‖A0‖(1-α)/ατα‖eA0τ‖1/α+∑l=1N1|Γ((1-l)α)|‖A0‖-lτ(1-l)α-1+O(‖A0‖-(N+1)τ-Nα-1)∑l=1N1|Γ((1-l)α)|)‖Ai‖‖Ψα(t-τ-hi)‖dτ
for all j∈k-1¯∪{0} as t→∞ if 0<α<2, for any N∈Z+, and(5.16)‖Ψαj0(t)‖≤1α(‖A0‖1/α)-j‖eA0t‖1/α+∑l=1N1|Γ(j+1-αl)|‖A0‖-ltj-lα+O(‖A0‖-(N+1)tj-(N+1)α)+∑i=1p∫0t(1α‖A0‖(1-α)/α‖eA0τ‖1/α+∑l=1N1|Γ((1-l)α)|‖A0‖-lτ(1-l)α-1+O(‖A0‖-(N+1)τ-Nα-1)∑l=1N1|Γ((1-l)α)|)‖Ai‖‖Ψαj0(t-τ-hi)‖dτ,‖Ψα(t)‖≤1α‖A0‖(1-α)/α‖eA0t‖1/α+∑l=1N1|Γ((1-l)α)|‖A0‖-ltj-lα+O(‖A0‖-(N+1)tj-(N+1)α)+∑i=1p∫0t(1α‖A0‖(1-α)/α‖eA0τ‖1/α+∑l=1N1|Γ((1-l)α)|‖A0‖-lτ(1-l)α-1+O(‖A0‖-(N+1)τ-Nα-1)∑l=1N1|Γ((1-l)α)|)‖Ai‖‖Ψαj0(t-τ-hi)‖dτ
for all j∈k-1¯∪{0} as t→∞ if α≥2, for any N∈Z+.For further discussion, note that there exists a set of linearly independent continuously differential real functions{αi:R0+→R,i∈ν-1¯∪{0}}, where ν is the degree of the minimal polynomial of  any square real matrix A0 such that:(5.17)eA0t=∑i=0ν-1αi(t)A0i=∑j=0ν∑i=0νi-1kijtieλjt;∀t∈R0+
(see, e.g., [4, 5]), where kij∈R; i∈ν-1¯∪{0}, j∈ν¯∪{0}, σ(M)∶={λi∈C:det(λiIn-A0)=0} is the spectrum of A0 defined by the  set of eigenvalues λiof M of respective index νi (i.e., the  multiplicity of λiin the minimal polynomial of A0) and algebraic multiplicity μi (i.e., the multiplicity of λi in the characteristic polynomial of A0) so that  n=∑i=1nni≥ν=∑i=1nνi with n being the order of A0 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) Ifk=α=1 (the particular standard bon-fractional case) then (3.1) is globally Lyapunov stable independent of the delays if
(5.18)‖(1β1A1,1β2A2,…,1βpAp)‖2≤-μ2(A0)
requiring for the ℓ2-matrix measure of A0 to fulfil μ2(A0)∶=(1/2)λmax(A0+A0T)≤0, for some βi∈R+(i∈p̅) subject to ∑i=1pβi2=1, [6]. Also,
(5.19)Ψ100(t)=eA0t-A0-1t-1+O(A0-1t-1)+∑i=1p∫0t(eA0τ-A0-1τ-1+O(A0-1τ-1))AiΨ100(t-τ-hi)dτast⟶∞
is bounded provided that A0 is non-singular with eA0t being of the form  (5.17) if (5.18) holds and  then the  unforced solution:
(5.20)xα(t)=∑j=0k-1(Ψαj0(t)xj0+∑i=1p∫0hiΨαj0(t-τ)φj(τ-hi)dτ).
Is bounded for all time. Furthermore,
(5.21)‖Ψ100(t)‖≤‖eA0t‖+‖A0‖-1t-1+O(‖A0‖-2t-2)+∑i=1p∫0t(‖eA0τ‖+‖A0‖-1τ-1+O(‖A0‖-1τ-1))×‖Ai‖‖Ψ100(t-τ-hi)‖dτast⟶∞
if (5.18) holds irrespective of A0 being singular or non-singular. If, in addition, μ2(A0)<0 and (5.18) holds with strict inequality then (3.1) is globally asymptotically Lyapunov stable independent of the delays and
(5.22)Ψ100(t)⟶∑i=1p∫0t(eA0τ-A0-1τ-1+O(A0-1τ-1))AiΨ10(t-τ-hi)dτ⟶0ast⟶∞.
(ii) Ifk=1 and α∈(0,1] the inequality (5.18) is strict then (3.1) is globally Lyapunov stable independent of the delays if μ2(A01/α)≤0 and
(5.23)‖(1β1A1,1β2A2,…,1βpAp)‖2≤-μ2(A01/α)
provided that A0 is non-singular and A01/α exists. Also, then (3.1) is globally asymptotically Lyapunov stable independent of the delays if, in addition, μ2(A01/α)<0 and
(5.24)‖(1β1A1,1β2A2,…,1βpAp)‖2<|μ2(A01/α)|,(5.25)Ψα00(t)=1αe(A01/α)t-1Γ(1-α)A0-1t-α+O(A0-2t-2α)+∑i=1p∫0t(1αA0(1-α)/αe(A01/α)τ-A0-1τ-1+O(A0-(N+1)τ-α-1))×AiΨα00(t-τ-hi)dτ⟶0ast⟶∞.
If either A0 is singular or A01/α does not exists then (5.25) is replaced by a corresponding less than or equal to  relation of norms with the replacements A0→∥A0∥, A0-1→∥A0∥-1, and eA0t→∥eA0t∥.
(iii) Assume thatJA0=JA0d+J̃A0 is the canonical real form of A0 (in particular, its  Jordan form if all the eigenvalues are real)  with JA0d being diagonal and J̃A0being off diagonal such that the above decomposition is unique with A0=T-1JA0T where T is a unique non-singular transformation matrix. Then, the Caputo fractional differential system (3.1) is globally Lyapunov stable  independently of A01/α to exist or not by replacing μ2(A01/α)→μ2(JA0d1/α) in (5.23) by
(5.26)‖(1β0T-1J̃A0T,1β1T-1A1T,1β2T-1A2T,…,1βpT-1ApT)‖2≤|μ2(JA0d)|1/α
with μ2(JA0d1/α)≤0 for some set of numbers βi∈R+(i∈p̅∪{0}) satisfying  ∑i=0pβi2=1. The fractional system is globally asymptotically Lyapunov stable for one such a set of real  numbers if μ2(JA0d1/α)<0, what implies that |argλ|<απ/2, for all λ∈σ(A0),  and
(5.27)‖(1β0T-1J̃A0T,1β1T-1A1T,1β2T-1A2T,…,1βpT-1ApT)‖2<|μ2(JA0d)|1/α.Proof.
It turns out thatxα(t) is bounded for all time so that (3.1) is globally Lyapunov stable if ∥Ψαj0(t)∥ is bounded; forallj∈k-1¯∪{0} for all t∈R0+ for any bounded functions of initial conditions φj:[-h,0]→Rn; for all j∈k-1¯∪{0} with φj(0)=xj(0)=xj0.  If, in addition, ∥Ψαj0(t)∥→0 as t→∞ then xα(t)→0 as t→∞ so that (3.1) is globally asymptotically Lyapunov stable and the solution (5.20) is bounded for all time. Thus, if k=α=1 (the particular standard bon-fractional case) then (3.1) is globally Lyapunov stable if ∥Ψ100(t)∥ is bounded for all t∈R0+. A sufficient condition independent of the delays is that (5.18) holds requiring trivially for the ℓ2-matrix measure of A0 to fulfil μ2(A0)∶=(1/2)λmax(A0+A0T)≤0, where the for some βi∈R+(i∈p̅) subject to ∑i=1pβi2=1, [27].  Equation (5.19) follows from (5.11)  after inspection for N=1 and it is bounded  as t→∞ and since otherwise the  global stability property (5.18) would fail contradicting its sufficient condition for j+1=k=α=1. Equation (5.20) follows from (5.14) for j+1=k=α=N=1 irrespective of A0 being singular or non-singular and of the fact of A01/α to exist or not. Equation (5.21) follows from (5.19) since μ2(A0)<0 implies that A0 is a stability matrix then Re(λ)<0; forallλ∈σ(A0)  and,  furthermore, Ψ100(t)→0, and the unforced solution xα(t)→0, as t→∞ from the strict inequality guaranteeing global asymptotic stability independent of the delays, namely, ∥(1/β1)A1,(1/β2)A2,…,(1/βp)Ap∥2<|μ2(A0)|. Property (i) has been proven. Property (ii) has a similar proof for α∈(0,1], k=1 by replacing A0→A01/α. Property (iii) follows by using the matrix similarity transformation A0=T-1JA0T=T-1(JA0d+J̃A0)T and using the homogeneous transformed Caputo fractional differential system from (3.1):
(5.28)(D0+αzC)(t)=(D0+αTxC)(t)=∑i=0pAiTx(t-hi)⟺(D0+αxC)(t)=∑i=0pT-1AiTx(t-hi)=T-1A0Tx(t)+∑i=1pT-1AiTx(t-hi)=T-1JA0dTx(t)+∑i=0pT-1A̅iTx(t-hi),
where z(t)=Tx(t); for all t∈R0+, h0=0 plays the role of an additional delay. A̅0=J̃A0 and A̅i=Ai(i∈p̅) by noting also that since (JA0d+JA0d*) is diagonal with real eigenvalues by construction, one  has
(5.29)|μ2(JA0d1/α)|=|12λmax⁡(JA0d1/α+JA0d1/α*)|=|λmax⁡(JA0d1/α)|=|Reλmax⁡(JA0d1/α)|=|Reλmax⁡1/α(JA0d)|=|Reλmax⁡1/α(A0d)|=|μ2(JA0d)|1/α.
Then, the proof is similar to that of the related part of Property (ii). Note also that μ2(JA0d1/α)<0, implies that
(5.30)Reλ1/α<0⟺arg(λα)∈(-π2,π2);∀λ∈σ(A0)⟺|argλ|>απ2,∀λ∈σ(A0).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 A0→T-1JA0dT and Ai→T-1A̅iT with A̅i(i∈p̅∪{0}) defined in the proof of Theorem 5.1(iii). Theorem 5.1 establishes  that for any stability matrix A0, the asymptotic stability condition of sufficient type  is as follows:(5.31)‖(1β0T-1J̃A0T,1β1T-1A1T,1β2T-1A2T,…,1βpT-1ApT)‖2<|μ2(JA0d)|1/α
provided that μ2(JA0d1/α0)<0 extends from α=α0≤1, (in particular, from  the standard nonfractional differential system α=α0=1) to any α∈(0,α0] provided that arg(λ/α)∈(-π/2,π/2) in the clockwise sense, or equivalently, if |argλ|>απ/2, for all λ∈σ(A0), since
(5.32)‖(1β0T-1J̃A0T,1β1T-1A1T,1β2T-1A2T,…,1βpT-1ApT)‖2<|μ2(JA0d)|1/α0≤|μ2(JA0d)|1/α,∀α∈(0,α0].
Note that the global Lyapunov’s stability conditions (5.23) and (5.26) with nonpositive measures μ2(JA0d1/α) 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 μ2(JA0d1/α), 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 φjis required to be identically zero on its definition domain for k-1¯∪{0}∋j<α-1(α≥2) in order that the Γ-functions be positive (note that Γ(x) is discontinuous at zero with an asymptote to -∞ as x→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 μ2(JA0d1/α). Assume also that φj:[-h,0]→Rn are any admissible functions of initial conditions for k-1¯∪{0}∋j≥α-1 while they are identically zero if  k-1¯∪{0}∋j<α-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].



---

*Source: 161246-2011-03-23.xml*