Asymptotic Stability of the Solutions of Neutral Linear Fractional System with Nonlinear Perturbation

: In this article existence and uniqueness of the solutions of the initial problem for neutral nonlinear differential system with incommensurate order fractional derivatives in Caputo sense and with piecewise continuous initial function is proved. A formula for integral presentation of the general solution of a linear autonomous neutral system with several delays is established and used for the study of the stability properties of a neutral autonomous nonlinear perturbed linear fractional differential system. Natural sufﬁcient conditions are found to ensure that from global asymptotic stability of the zero solution of the linear part of a nonlinearly perturbed system it follows global asymptotic stability of the zero solution of the whole nonlinearly perturbed system.


Introduction
Currently, a variety of scientific fields are successfully using the latest advances in fractional calculus and fractional differential equations. For a good introduction in the theory of fractional calculus and fractional differential equations see Kilbas et al. [1], Kiryakova [2] and Podlubny [3]. The distributed order fractional differential equations is discussed in Jiao et al. [4] and for an application oriented exposition see Diethelm [5]. We refer also the monograph of Stamova, Stamov [6] where impulsive fractional differential and functional differential equations as well as several applications are considered.
It is well known that the stability of a process is the ability of the process to withstand previously unknown, small influences (perturbations). If such perturbations do not substantially change the process, then it is called stable. We emphasize that this property proves to be extremely important and becomes an "evergreen" research topic. As in the integer case, the study of the stability of fractional differential equations and systems with delay is more complicated compared with fractional differential equations and system without delay. We point out that this is due to the fact that, in fractional delay differential equations, the dependence on the past evolution history of the processes described by such equations is inspired by two sources. First of them is the impact conditioned by the delays and the other one the impact conditioned from the availability of Volterra type integral in the definitions of the fractional derivatives, i.e., the memory of the fractional derivative. It must be noted that the first of them (conditioned by the delays) is independent from the derivative type (integer or fractional). Different types fractional differential equations and systems with delays (retarded and neutral) or without delays are studied for several types of stability. As works related to this theme we refer to [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].
In this article, first a general case of nonlinear delayed fractional system with linear neutral part and variable delays is considered. The fractional derivatives of the system are in Caputo sense with incommensurate orders α k ∈ (0, 1), k = 1, . . . , n. The incommensurate order of the fractional derivatives means that, unlike many fractional systems studied, the order of the fractional derivative is not the same for the whole system, and moreover, the different orders of the fractional derivatives are not rational numbers, which would allow a common denominator to be found (such approach has also been widely used in some studies). For this type systems, we prove existence and uniqueness of the solutions of an initial problem (IP) with piecewise continuous initial conditions. We know only a few results for Cauchy problem for fractional delay differential equations with initial functions which are not continuous (see [24][25][26]).
Then we have two main goals. First of them is to obtain sufficient conditions which guarantee that the zero solution of a neutral linear system with nonlinear perturbation is globally asymptotically stable if the zero solution of the unperturbed neutral linear system is globally asymptotically stable. The second one is to study the influence of the memory on the asymptotic nature of the solutions of the these systems, which is generated by the fractional derivatives and the time delays in the systems. Since the conditions and the obtained results are similar as these in the case of delayed systems with integer derivatives we can conclude that the influence from the memory generated by the time delays in the systems has more determining influence for the evolution of the process in compare with this generated by the fractional derivatives.
It must be noted that for the study of the stability properties described above, a formula for integral representation of the general solution of a linear autonomous neutral system with several delays is proved. For papers, related to such problems we refer to [24,[27][28][29].
The paper is organized as follows. In Section 2 we give definitions and needed properties of Riemann-Liouville and Caputo fractional derivatives and introduce some notations. In Section 3 we prove existence and uniqueness of the solutions of the initial problem for neutral nonlinear differential system with incommensurate order Caputo fractional derivatives and with piecewise continuous initial function. In Section 4 we establish a formula for integral presentation of the general solution of a linear autonomous neutral system with several delays which is needed in our investigations below. Note that the obtained result are an immediate generalization of the results obtained in [27]. Section 5 is devoted to the study of a neutral autonomous nonlinear perturbed linear fractional differential system in the case of Caputo type derivatives with incommensurate differential orders. Using the formula obtained in the previous section, some natural sufficient conditions are found to ensure that from global asymptotic stability of the zero solution of the linear part of a nonlinearly perturbed system it follows global asymptotic stability of the zero solution of the whole nonlinearly perturbed system.

Preliminaries
Let γ ∈ (0, 1) be an arbitrary number and denote by L loc 1 (R, R) the linear space of all locally Lebesgue integrable functions f : R → R. Then for a ∈ R, each t > a and f ∈ L loc 1 (R, R) the left-sided fractional integral operator, the left-sided Riemann-Liouville and Caputo fractional derivative of order γ are defined by respectively (see [1]). We will use the following relations (see again [1]): Concerning the Laplace transform (L f )(p) = ∞ 0 e −pt f (t)dt, p ∈ C we need the properties: In this article we will use only one-side Laplace transform. The main criterion that we use for the existence of a Laplace transform is the exponential boundedness of the functions. For more details on Laplace transform see [30].
Everywhere below we will use the notations R + = (0, ∞), . . , n}, n ∈ N, n 0 = n ∪ {0}, I, Θ ∈ R n×n denote the identity and zero matrix respectively and 0 ∈ R n is the zero element. For We will use also the notations As usual for arbitrary fixed h > 0 a vector function

Existence and Uniqueness of the Solutions of the Cauchy Problem for Neutral Nonlinear Fractional Differential System
Consider the nonlinear delayed system of neutral type with incommensurate Caputo fractional derivatives or described in more detailed form . , x n t (θ)), x k t (θ) := x k (t + θ) for t ∈ J a , −h ≤ θ ≤ 0, h > 0 be an arbitrary fixed constant. As previously explained we consider x k t (θ) for every fixed t ∈ J a as the restriction of the function x k (t) on the interval [t − h, t] (see [31,32]).
Introduce for arbitrary Φ ∈ C * the following initial condition for both types of delays i.e., for each k ∈ n we have that x k a (θ) = x k (a + θ) = φ k (θ) for −h ≤ θ ≤ 0 and For the neutral part of the system (1) we say that the conditions (A) are fulfilled if the following conditions hold: (A1) The matrices A l (t) = {a l kj (t)} n k,j=1 ∈ C(J a , R n×n ) for every l ∈ r . (A2) The delays τ l (t) ∈ C(J a , R + ), τ l (a) > 0 and sup t∈J a τ l (t) ≤ h for every l ∈ r .
(A3) The set S * Φ = {t ∈ J a |t − τ l (t) ∈ S Φ , l ∈ r } do not have limit points. Consider the following auxiliary system or described in more detailed form for k ∈ n , is a solution of the IP (1) and (2) or of the IP (3) and (2) in [a, a + M] (J a ), if it satisfies the system (1) or respectively (3) for all t ∈ (a, M] (t ∈ (a, ∞)) and the initial condition (2) too.
We say that for the vector valued functional F : E * → R n the ((H)/Caratheodory/conditions are fulfilled in E * if the following conditions hold: (H1) For almost all fixed t ∈ J a the function (t, Ψ) → F(t, Ψ) is continuous in arbitrary Ψ ∈ C * and for each fixed function Ψ ∈ C * the function (t, Ψ) → F(t, Ψ) is Lebesgue measurable and locally bounded for t ∈ J a .

Remark 1.
Note that the Lipschitz conditions (4) in (H2) imply that for each t ∈ J a we have F(t, 0 T ) ≡ 0. Furthermore, the function ∈ L loc 1 (J a , R + ) in (H2) can depend from the neighborhood of the chosen point (t, Ψ) ∈ E * . For more details about Lipschitz functions see [33]. Proof. The proof is almost the same as the proof of the Lemma 1 in [34] for the case of continuous initial function but for completeness we will sketch it.
Let X(t) = (x 1 (t), . . . , x n (t)) T be a solution of the IP (1) and (2) in J a . Then condition (H1) implies that F(t, X T t ) is Lebesgue integrable function. Applying the operator D −α k a+ , k ∈ n to both sides of (1) and using formula (c) we obtain that for the left side of (1) the following equality holds where the constant c k Φ is calculated by the use of the initial conditions (2). Then from (1) and (2) and (5) it follows that X(t) is a solution of the IP (3) and (2).
Conversely if X(t) is a solution of the IP (3) and (2) then we apply the operator D α k a+ , k ∈ n , to both sides of (3) and taking into account (b) and (5) we obtain that X(t) is a solution of the IP (1) and (2).
For arbitrary fixed Φ ∈ C * we introduce the following set and for arbitrary M ∈ R + the sets a be arbitrary and introduce in E * a the following distance function It is simply to check that the sets E * M and E * a are complete metric spaces in respect to the introduced distance function. Theorem 1. Let the following conditions be fulfilled:

1.
For the vector valued functional F : E * → R n the conditions (H) hold in E * and the conditions (A) hold too.

2.
The initial function Φ ∈ C * has at most one jump point t Φ ∈ [a − h, a] and is right continuous on S Φ .
Then there exists M 0 ∈ R + such that the IP (3) and (2) has a unique solution in the interval [a, a Condition (A2) implies that there exists c * ∈ R + such that for t ∈ [a, a + c * ] and all l ∈ r the inequalities t − τ l (t) < a hold. Without loss of generality we can assume that c * ≤ min{1, h}. Let M ∈ (0, c * ], t * ∈ [a, a + M] be arbitrary and then for an arbitrary function W t * = (w 1 t * , . . . , w n t * ) T ∈ J * M define the operator (RW t * ) = ((Rw 1 t * ), . . . , (Rw n t * )) T point wise for every t ∈ [t * − h, t * ] as follows: or for k ∈ n in more detailed form: For the second addend in (6) from Condition 1 of the theorem it follows that for each t ∈ (a, t * ]. Then from Condition 1 of the theorem, (8) and (9) it follows that the second addend in the right side of (6) is a continuous function for t ∈ (a, t * ] and hence (6) implies that the function From (6), (7), (8) and (10) for every t ∈ [a, a + M] we obtain that . . , w n s )|ds. (s) and then from (11) where , α m = min(α 1 , . . . , α n ) and α M = max(α 1 , . . . , α n ).
Then as in the former case (a) we can prove that there exists M 0 ∈ (0, c * ] such that the operator R defined by (6) Since t − τ l (t) are continuous functions at a and t Φ < a for all l ∈ r with l / ∈ {l j 1 , . . . , l j p } we can conclude that there exists c * ∈ (0, c * 1 ) such that for t ∈ [a, a + c * ] the inequality min t∈[a,a+c * ] (t − τ l (t)) > t Φ holds. Then the same way as in the proof of point (a) above, we can obtain that there exists M 0 ∈ (0, c * ] such that the operator R defined by (6)-(8) is contractive in E * M 0 .

Remark 2.
Note that from Theorem 1 it follows that any solution of the IP (3) and (2) is unique on the interval where this solution does exist. That's mean if there exist two solutions X 1 (t), X 2 (t) of the IP (3) and (2) with intervals of existence [a, a + M 1 ] and [a, the solution X 2 (t) is a continuation of X 1 (t).

Remark 3.
It is not hard to check that the proof of Theorem 1 remains useful in the essential more general case with finitely many first kind jumps of the initial function Φ ∈ C * when the intersection S * Φ ∩ T a = ∅ holds.
The aim of the next corollary is to study the important case of the, I. when the right end of the initial interval does not coincide with the lower terminal of the fractional derivatives.
Let X M 0 (t) be the unique solution of IP (3) and (2) in the interval [a, a + M 0 ]. Consider the initial condition for the system (3) with shifted initial point t 0 = a + M 0 and initial function X M 0 (t), t ∈ [a − h, a + M 0 ] as follows: , is a solution of the IP (1) and (13) or of the IP (3) and (13) in [t 0 , t 0 + M] (J t 0 ), if it satisfies the system (1) or respectively (3) for all t ∈ (t 0 , M] (t ∈ (t 0 , ∞)) and the initial condition (13) too.

Remark 4.
Let X M t 0 (t) be the unique solution of IP (3) and (2) in the interval [a, a + M 0 ]. Then if we choose t 0 = a + M 0 as initial point and take X M t 0 (t) as initial function in the interval [a − h, t 0 ] for the IP (3) and (13), then using the solution of IP (3) and (13) (if there exists) we can define a prolongation of X M t 0 (t) as solution of the IP (3) and (2).
Note that the most complicated case is when t 0 < a + h and t Φ = a. Below we will consider only this case. Corollary 1. Let the following conditions hold.

1.
The conditions of Theorem 1 hold.
Then there exists M 1 > 0 such that the IP (3) and (13) has a unique continuous solution in the interval Proof. The proof is almost the same as the proof of Theorem 1 but for completeness we will sketch it.
As above for arbitrary fixed Φ ∈ C * we introduce the following set and for arbitrary M > 0 the sets a we see that G ∈ C * for every fixed t ∈ J a . Then for arbitrary M > 0 we have that Let M > 0 and t * ∈ [t 0 , t 0 + M] be arbitrary and then for every function . . , (R n w n t * )) T for t ∈ (t * − h, t * ] as follows: Define the operator R for t ∈ (t 0 , t * ) with (6); for t = t * with (8) and Note that (14) is similar condition as (7) but with other initial point and initial function.
Consider the set T t = {t − τ l (t) | l ∈ r } for t ∈ [a, a + h].
Consider also the set T t 0 = {t 0 − τ l (t 0 ) | l ∈ r } and let a / ∈ T t 0 , i.e., S * Φ ∩ T t 0 = ∅. Then as in the case (b) of Theorem 1 from conditions (A) it follows that there exists c * ∈ R + , such that for t ∈ [t 0 , t 0 + c * ] we have S * Φ ∩ T t = ∅, i.e., a / ∈ T t . Thus for t ∈ [t 0 , t 0 + c * ] we have that t − τ l (t) is a continuous function for each l ∈ r .
Let M ∈ (0, c * ] be arbitrary. Then for every t * ∈ [t 0 , t 0 + M] and each W t * ∈ J * M from (14) and from Condition 1 of Theorem 1 for t ∈ [t * − h, t 0 ] we have that (15) and hence the second addend in (6) is a continuous function for t ∈ [a, t 0 ] (right continuous at a). Moreover, from Condition 1 of Theorem 1 and (8) it follows that the second addend in the right side of (6) is a continuous function for t ∈ [t 0 , t * ] and thus (6) implies that the function (RW t * )(t) is continuous for t ∈ [t 0 , t * ] too. Since from (15) and (6) it follows that X M 0 (t) is continuous at t 0 then we can conclude that R(J * Then the same way as in the proof of Theorem 1 we obtain Consider the case when a ∈ T t 0 . Then S * Φ ∩ T t 0 = ∅. Then as in the case (c) of Theorem 1 from conditions (A) it follows that there exist some numbers l j 1 , . . . , l j p , 1 ≤ p ≤ r, such that t 0 − τ l j i (t 0 ) = a, 1 ≤ i ≤ p. For every ε ∈ (0, δ Φ ), where δ Φ is the same as in Theorem 1, since Φ(t) is right continuous at a then there exists c * 1 ∈ R + , such that for t ∈ [t 0 , t 0 + c * 1 ] we have |Φ(a) − Φ(t − τ l j i (t))| < ε. Thus for t ∈ [t 0 , t 0 + c * 1 ] we have that t − τ l j i (t) ≥ a. Since for all l ∈ r with l / ∈ {l j 1 , . . . , l j p } the functions t − τ l (t) are continuous at t 0 with t 0 − τ l (t 0 ) = a, then we can conclude that there exists c * ∈ (0, c * 1 ) such that for l ∈ r with l / ∈ {l j 1 , . . . , l j p } we have that S * Φ ∩ {t − τ l (t) | t ∈ [t 0 , t 0 + c * ]} = ∅. Thus for l ∈ r with l / ∈ {l j 1 , . . . , l j p } we have that a / ∈ {t − τ l (t) | t ∈ [t 0 , t 0 + c * ]} and hence the functions t − τ l (t) are continuous for these l and t ∈ [t 0 , t 0 + c * ]. Then the same way as in the proof of the former case above, we can obtain that there exists M 1 ∈ (0, c * ] such that the operator R defined by (6), (8) and (14) is contractive in E * M 1 .

Theorem 2. Let the conditions of Theorem 1 hold. Then the IP (3) and (2) has a unique solution in J a .
Proof. According Theorem 1 there exists M 0 > 0 such that the IP (3) and (2) has a unique solution in [a, a + M 0 ]. Denote by X max (t) = (x max 1 (t), . . . , x max n (t)) the maximal solution of the IP (3) and (2) and assume that the interval of existence J max is closed from right, i.e., J max = [a, a + M max ] and X max (t) is a continuation of every other solution of the IP (3) and (2). Then applying Corollary 1 with initial point M max and initial function X max (t) we obtain a prolongation of X max (t) which is a contradiction. Thus we conclude that the interval of existence has the form J max = [a, a + M max ).
Let we assume that M max < ∞. Then we have two cases : either a + M max − τ l (a + M max ) = a for every l ∈ r , or there exist some numbers l j 1 , . . . , l j p , 1 ≤ p ≤ r, such that a + M max − τ l j i (a + M max ) = a, 1 ≤ i ≤ p. Let consider the case when a + M max − τ l (a + M max ) = a for every l ∈ r . Then the right side of (3) is continuous and passing to limit in the both sides of (3) for t → a + M max − 0 we obtain that (3) holds for t = a + M max . Therefore we are obtained a solution which is a prolongation of X max (t) since it has as interval of existence [a, a + M max ] which is a contradiction and hence M max = ∞ in this case.
Let there exist some numbers l j 1 , . . . , l j p , 1 ≤ p ≤ r, such that a + M max − τ l j i (a + M max ) = a, 1 ≤ i ≤ p. Then we have that and since Φ ∈ C * , then the right side of (17) has finite limit. Therefore the right side of (3) can be prolonged as continuous function at t = a + M max as well as the left side and therefore (3) holds for t = a + M max too. Thus M max = ∞ in this case too.

Integral Representation of the Solution of the, I. for Autonomous Linear Neutral Fractional System
The aim of this section is to obtain an integral representation of the solutions of autonomous linear fractional neutral system with Caputo type derivatives and multiple delays introduced below (see (19)). The obtained representation will be essentially used in the next Section 5.
As usual a vector valued function X T (t) = (x 1 (t), . . . , x n (t)) ∈ C(R + , R n ) will be called exponentially bounded, if for t ∈ R + we have that |X(t)| ≤ Ce tγ for some C > 0 and γ ∈ R.
Consider an autonomous linear neutral fractional system with derivatives in Caputo sense and multiple delays in the following form and the homogeneous one where A i , B i ∈ R n×n , σ i ∈ (0, σ], σ 0 = 0, τ l ∈ (0, τ], σ, τ ∈ R + , i ∈ m , l ∈ r , X, F : R + → R n , h = max(σ, τ). Consider the following initial conditions for the systems (18) or (19): Let s ∈ R + be an arbitrary fixed number and consider the following matrix, I. for t ∈ J s D α where Q(t, s) = {γ kj (t, s)} n k,j=1 : [s, ∞) × R + → R n×n and initial condition Definition 3. For each s ∈ R + the matrix valued function t → Q(t, s) is called a solution of the IP (21), (22) for t ∈ J s = [s, ∞) if Q(·, s) is continuous in t on J s and satisfies the matrix Equation (21) for t ∈ (s, ∞), as well as the initial condition (22) too.
In the case when s = 0, the matrix Q(t) = Q(t, 0) will be called the fundamental (or Cauchy) matrix of a system (19).

Remark 5.
Note that from Theorem 2 it follows that the matrix IP (21) and (22) has a unique solution. Moreover, from Theorem 2 in [34] it follows that the IP (18) and (20) has a unique continuous solution for each Φ ∈ C = C([−h, 0], R n ) and locally bounded F ∈ L loc 1 (R + , R n ). It must be also noted that for the Equations (18) and (19) the conditions (A) are fulfilled.
The next results are an immediate generalization of the results obtained in [27] . Theorem 3. The fundamental matrix Q(t) of (19) is exponentially bounded and has the following representation where is the characteristic matrix of (19) (see [34]).
Proof. Let us assume that every column of the fundamental matrix Q(t) of (19) is exponentially bounded, i.e., is O(e δt ) in general for some δ > 0 . Then we can correct apply the Laplace transform to both sides of (21) and similar as in the proof of the corresponding result in [27] we obtain that the representation (23) holds. Hence the matrix Q(t) = Q(t, 0) is a solution of IP (21) and (22) for s = 0 . Since the IP (21) and (22) in virtue of Theorem 2 has a unique solution then we obtain that the matrix Q(t) defined by (23) is this unique solution. Since the real parts of the roots of the characteristic equation detG(p) = 0 are uniformly bounded from above, then from the representation (23) it follows immediately that the fundamental matrix Q(t) of (19) is exponentially bounded.

Theorem 4.
For every Φ ∈ C the corresponding unique solution X Φ (t) of the IP (19) and (20) can be represented in the following form: where Q(t) is the fundamental matrix of (19).
Proof. Since Theorem 3 implies that the fundamental matrix Q(t) of (19) is exponentially bounded, then from (24) it follows that X Φ (t) is exponentially bounded too. Substituting X Φ (t) in (19) and applying the Laplace transform to both sides of (19) we obtain that where Introduce the functions: for every i ∈ m and l ∈ r . Then using for each l ∈ r the substitution s = t + τ l we obtain The same way we obtain Taking into account (26)-(28) we receive and applying to both sides of (29) the inverse Laplace transform we have For every l ∈ r after simple calculation we obtain and from (31) for the second addend in (30) it follows that Since for the fourth and fifth addends in the right side of (30) we have that (33) and then substituting in (30) the results from (32) and (33) we obtain which completes the proof. Theorem 5. Let the function F ∈ L loc 1 (R + , R n ) be exponentially bounded. Then the solution X F (t) of the IP (18) and (20) with initial function Φ(t) ≡ 0, t ∈ [−h, 0] has the following representation: where Q(t) is the fundamental matrix of the system (21).
Proof. The proof of this result is almost the same as the proof of the corresponding result in [27] and will be omitted.

Corollary 2.
Let the function F ∈ L loc 1 (R + , R n ) be exponentially bounded. Then for every initial function Φ ∈ C the corresponding unique solution X F Φ (t) of the IP (18) and (20) has the following integral representation: where Q(t) is the fundamental matrix of system (19).
Proof. Let Φ ∈ C be an arbitrary initial function and let the function X Φ (t) be the unique solution of IP (19) and (20) with Φ ∈ C and let X F (t) be the unique solution of IP (18) and (20) with initial function Φ(t) ≡ 0, t ∈ [−h, 0] for arbitrary exponentially bounded function F ∈ L loc 1 (R + , R n ). Then according the superposition principle the function X F Φ (t) = X Φ (t) + X F (t) is the unique solution of IP (18) and (20).

Asymptotic Stability of a Nonlinear Perturbed Fractional System with Neutral Autonomous Linear Part
Consider the neutral nonlinear perturbed system i.e., where X T (t) = (x 1 (t), . . . , x n (t)) ∈ C(R + , R n ), W : E → R n , E = R + × C, W T = (w 1 , . . . , w n ), w k : E → R, k ∈ n and which neutral linear part coincides with the system (19).
For the system (36) introduce the following initial condition Remark 6. It is well known that the system (36) is a partial case of the system (1). Everywhere below we will assume that the initial point is a = 0. Theorem 6. Let the following conditions be fulfilled: 1.
The conditions (A) hold.

2.
For the vector valued functional W : E → R n in the right side of the perturbed system (36) the conditions (H) hold for each (t, Ψ) ⊂ E.
Then for every fixed initial function Φ ∈ C the IP (36) and (37) has a unique solution in R + .
Proof. The statement of Theorem 6 follows immediately from Theorem 2.

Definition 4.
We say that the vector valued functional W : E → R n is exponentially bounded in C(R + , R n ) if for every X T (t) ∈ C(R + , R n ) there exist constants C X ∈ R + , γ X ∈ R (i.e., the constants can depend from X) such that for the the function F(t) = W(t, X T t ) holds |F(t)| ≤ C X e γ X t for t ∈ R + .

Definition 5.
The zero solution of the system (18), (19) or (37) is said to be: (a) Stable (uniformly) iff for any ε > 0 there is a δ(ε) > 0 such that for every initial function Φ ∈ C with ||Φ|| < δ the corresponding solution X(t) satisfies for each t ∈ R + the inequality |X(t)| ≤ ε. (b) Locally asymptotically stable (LAS) iff there is a ∆ ⊂ C such that for every initial function Φ ∈ ∆, the relation lim t→∞ |X(t)| = 0 holds for the corresponding solution X(t).
(c) Globally asymptotically stable (GAS) iff for every initial function Φ ∈ C, for the corresponding solution X(t) we have that lim t→∞ |X(t)| = 0.
The next simple lemma plays an important role in the proof of the main result in this section.

Lemma 2.
Let Q(t) = Q(t, 0) is the fundamental (or Cauchy) matrix of system (19) in the case when s = 0 and the zero solution of the system (19) is globally asymptotically stable (GAS). Then for every β ∈ (0, 1) we have that lim Proof. In virtue of Theorem 3 the fundamental matrix Q(t) of (19) has the following representation Q(t) = (L −1 [I α−1 (p)G −1 (p)](t). Then since the zero solution of the system (19) is GAS it follows that lim t→∞ Q(t) = Θ and all eigenvalues of the characteristic matrix G(p) of (19) belong to C − . Applying to the matrix C D β 0+ Q(t) the Laplace transform we obtain (L C D β 0+ Q(t))(p) = I β (L Q(t))(p) − I β−1 and hence we have that the function p(I β (L Q(t))(p) − I β−1 ) = I 1+β (L Q(t))(p) − I β is an entire function for p ∈ C + . Then taking into account that lim p→0 (I 1+β (L Q(t))(p) − I β ) = Θ (note that for p ∈ C + the function (L Q(t))(p) = I α−1 (p)G −1 (p)](t) is bounded), we can apply the final value theorem and hence lim The aim of the next theorem is to prove that if the zero solution of the system (19) (i.e., the linear part of system (36)) is GAS, then every solution X(t) of the IP (36), (37) with initial function Φ ∈ C is GAS. Theorem 7. Let the following conditions be fulfilled:

2.
The vector valued functional W : E → R n is bounded in C(R + , R n ).

3.
The zero solution of the system (19) is GAS.
Then every solution X(t) of the IP (36) and (37) with initial function Φ ∈ C is GAS.
Proof. Let for arbitrary initial function Φ ∈ C, X(t) be the unique solution of the IP (36) and (37). Substituting X(t) in (36) we obtain that where F(t) = W(t, X T t ) and hence according to Condition 2 of the theorem we have that |F(t)| ≤ C X and F(t) is piecewise continuous for t ∈ R + . Then from (38) and Corollary 2 we obtain that for X(t) the integral representation (35) holds, where Q(t) is the fundamental matrix of system (19). Under the conditions of the theorem we can apply the Laplace transform correct to both sides of (35) and after multiplying both sides of the received equality with p ∈ C + we obtain that p(L X(t))(p) = p(L Q(t))(p)C Φ + p(L D (p)(L Φ * τ l (t − τ l ))(p) It is clear that the right side of (39) is an entire function for p ∈ C + . Lemma 2 implies that the functions p(L Q(t))(p), p(L D 1 2 0+ Q(t))(p) and p(L D 1−α 0+ Q(t))(p) tends to 0 ∈ R n when p → 0 with p ∈ C + . Since the functions Φ * τ l , Φ τ l and Φ σ i , l ∈ r , i ∈ m 0 are piecewise continuous and bounded for t ∈ R, then we can conclude that the first five addends in the right side of (39) tend to 0 ∈ R n when p → 0 with p ∈ C + . From Condition 2 of the theorem it follows that F(t) is at least piecewise continuous for t ∈ R + and then Lemma 2 implies that the sixth addend tends to 0 ∈ R n when p → 0 with p ∈ C + too.
For the last addend we have that p(L D −α 0+ F(t))(p) = pp −α (L F(t))(p) = p 1−α (L F(t))(p) and hence the right side of the equality tends to 0 ∈ R n when p → 0 with p ∈ C + . Thus the right side of (39) tends to 0 ∈ R n when p → 0 with p ∈ C + . Then for p ∈ C + in virtue of the final value theorem we have that lim t→∞ X(t) = lim t→∞ p(L X)(p) = 0 ∈ R n .

Conclusions
In this article first the existence and uniqueness of the solutions of the initial problem in a general case for neutral nonlinear differential system with incommensurate order Caputo fractional derivatives and with piecewise continuous initial function is proved.
Then, knowing that such a solution exists, we look at the linear autonomous case and establish a formula for integral presentation of the general solution of a linear autonomous neutral system with several delays, which is an immediate generalization of the formula obtained in [27].
By the use of the obtained integral presentation of the general solution is studied a neutral autonomous nonlinearly perturbed linear fractional differential system in the case of Caputo type derivatives with incommensurate differential orders. Some natural sufficient conditions are found to ensure that from global asymptotic stability of the zero solution of the linear part of a nonlinearly perturbed system it follows global asymptotic stability of the zero solution of the whole nonlinearly perturbed system.
We hope that the results obtained will be useful both for future research and generalizations from a mathematical point of view, as well as for modeling of real-world phenomena.
Author Contributions: Conceptualization, H.K. and A.Z. Writing-Review and Editing, H.K. and A.Z. All authors contribution in the article are equal. All authors have read and agreed to the published version of the manuscript.