Backstepping sliding mode control for coupled delayed time fractional reaction diﬀusion systems subject to boundary input disturbances

In this paper we develop backstepping-based boundary sliding mode controllers for coupled time fractional reaction diﬀusion (FRD) systems governed by time fractional partial diﬀerential equations (PDEs) with time varying delay and input disturbances. Here, the diﬀusion coeﬃcients are spatially varying (nonconstant) and can be same or diﬀerent. To asymptotically stabilize the system, we ﬁrst use a backstepping transformation to convert the original dynamics into a target dynamics with a new manipulable input and the perturbation. Then, we employ the sliding mode algorithm to design this discontinuous controller to suppress disturbances. In this case, we obtain the combined back-stepping/sliding mode controller of the original dynamics, with which the closed-loop asymptotic stability is proved by the fractional Halanay’s inequality. Fractional numerical examples are given to verify the eﬃciency of the proposed synthesis.


Introduction
Fractional order models are performed to be more realistic and informative than the integer order ones. It is mainly due to the fact fractional derivatives [30] can be viewed as a powerful tool to describe memory and hereditary effects for modelling dynamical processes in a practical matter where integer order derivatives appear to fail. For the time fractional reaction diffusion (FRD) systems (anomalous reaction subdiffusion systems) modelled by time fractional partial differential equations (PDEs), it is obtained by replacing the firstorder on time t of reaction diffusion system with a fractional derivative of order α ∈ (0, 1]. This class system is not only full of theoretical value, but also has a lot of important applications. Examples include transport dynamics in systems governed by anomalous diffusion [23], combustion processes [36], biological systems [33], etc. Especially, fractional diffusion models has been used for modeling of the diffusion of COVID-19 [13,14], and the optimal vaccination and treatment policies of time fractional diffusion susceptible-infected-recovered (SIR) epidemic systems [17] could play a role in COVID-19. Additionally, the control problem of time FRD systems has been reported recently (see e.g. [16,31]). A first boundary control problem of a time FRD system with constant coefficients is developed in [18] using the backstepping approach [24]. The essence of this method in [24] is to utilize an invertible integral transformation to convert an original system into a target system. Subsequently, a fair amount of literatures (see e.g. [8,10]) on the counterparts of the case of spatially varying coefficients have emerged.
In the real-world problems, almost all control systems have coupling which makes the system complicated and needs to be addressed properly to avoid sys-tems unstable. Due to its practical significant importance, the boundary control/stabilization problems for coupled time FRD systems have also attracted much attention (see e.g. [10]). So far, these literatures have dealt almost only with the boundary disturbance-free control problems without time delay. But in practice, it is unavoidable to exist disturbances and time delay. For time FRD systems subject to input disturbance, the Mittag-Leffler stabilization problem has been considered in [37] by using the active disturbance rejection control (ADR-C) approach to design the boundary controller. Recent advances in the research area of time FRD systems with delay mainly focus on the periodicity and stability criteria [5], ADRC for stabilization problem [4], the existence of solutions [38], numerical methods for solutions [27], etc. However, to the best of our knowledge, the boundary control of coupled nonconstant-diffusivity time FRD systems with input matched disturbances and time varying delay has not been addressed yet.
Motivated by the aforementioned analysis, this paper aims to develop a backstepping-based sliding mode control strategy able to solve the asymptotic stability of a coupled nonconstant-diffusivity time FRD system with boundary input disturbances and time varying delay. Difficulties here lie in how to deal with the coupling, external disturbances, and the time varying delay. In this sense some special techniques (e.g. Filippov regularization approach [15], fractional Halanays inequality [20], sliding mode [7,34]) are employed.
The contribution is two-fold. First, we extend the combined backstepping/sliding mode boundary control method to the coupled nonconstant-diffusivity time FRD system subject to input disturbances and time varying delay. The resulting asymptotic stability is obtained by adopting the fractional Halanays inequality [20] besides using the fractional Lypunov method [28]. Then, we present the rigorous analysis of the well-posedness of the closed-loop system. Different from non-delayed cases in [37] and boundary disturbance-free cases in [9,10], herein we cope with the time varying delay with the help of the fractional Halanays inequality, and address the boundary disturbance by sliding mode control. Although the ADRC method adopted in [37] can be used to design the controller to compensate the input disturbance, it is restricted to observability of the extend state observer, and may be hard to adjust the corresponding parameters when we use it. In this sense, the sliding mode control method is an alternative to counteract the external disturbances of time FRD systems and this method has many advantages, such as robustness against disturbances and quick convergence properties.
This paper is organized as follows. In section 2, the problem illustration and preliminaries are provid-ed. Section 3 gives the boundary sliding mode control algorithm for the case of same diffusivity, including the well-posedness of the original system and the proof of the asymptotic stability. The results are generalized to the case of distinct diffusivity in Section 4. In Section 5 the numerical simulation is studied. Finally, Section 6 shows concluding remarks.
Notations: L 2 (0, 1) represents Hilbert space of a square integrable function u(x) with the norm ||u(·)|| = ( and λ min (·) stand for the largest and smallest eigenvalues respectively. In this work, we use f .., n, and I n×n is an n-dimensional identity matrix.

Problem illustration and preliminaries
Consider the below delayed system of coupled time fractional PDEs with external disturbances where x ∈ (0, 1), A(x) = diag(a 1 (x), ..., a n (x)) ∈ R n×n , n > 1, a i (x) denote the space-dependent diffusion coefficients that can be same or distinct, a max ≥ a i (x) ≥ a min > 0 with the positive constants a max , a min , i = 1, .., n, x ∈ [0, 1], U (x, t) ∈ R n stands for the state, represent the control input and the uncertain sufficiently smooth disturbance respectively. Here F (x, t) ∈ R n is an initial condition (IC), C 0 D α t U (x, t) denotes the Caputo fractional derivative of order α of U (x, t) with respect to t. Suppose that time varying delay τ (t) is continuously differential and satisfies In what follows, the system is actually coupled by the reaction coefficient Ψ (x) in the same diffusivity and the distinct diffusivity cases. For the system parameters, the below assumption is needed. Assumption 1. The coefficients of the system are sufficiently regular, in particular a i (x) ∈ C 2 [0, 1], ψ ij (x) ∈ C 1 [0, 1], i, j = 1, 2, ..., n.
Assumption 2. The initial function U (x, 0) in the IC (4) is compatible to the below BCs where the entries φ i (t) of Φ(t) are continuously differentiable and satisfy The open-loop system (1)-(4) (with V c (t) = 0) is potentially unstable, which depends on the values of the system coefficients and the boundary disturbance (see simulations for an illustrating example). The aim here is to design the boundary sliding mode controller by the backstepping method to guarantee asymptotic stability of the controlled system trajectories in the space [L 2 (0, 1)] n . In other words, the first task is to use the backstepping transformation to convert the original system into the target system with a new actuation and the perturbation. The second task is to employ the sliding mode algorithm to design this new control implement, which is used to reject the input perturbation. In this case, we obtain the asymptotic stabilization of the target dynamics, and therefore asymptotic stabilization of the original one.
The following definitions, lemmas and inequalities will be useful. [30]) The Caputo fractional derivative of order α is described by
Consider the following fractional differential inclusions where α ∈ (0, 1], T > 0, A 0 is the infinitesimal generator of an analytic semigroup of uniformly bounded linear operators on a Banach space H, X 1 (τ (t)) = X(t − τ (t)), and f is a multivalued map.
Definition 3. [1,2] Suppose that A 0 is a closed linear operator whose domain D(A 0 ) is defined on a Banach space H, and α > 0. Then we have that A 0 is the generator of an α-resolvent family if there exist ω ≥ 0 and a strong continuous function S α such that {λ α : In this sense S α (t) is named the α-resolvent family which is generated by A 0 .
Proof. From [32,Theorem 12.31] and Definition 3, we know that A 0 is the infinitesimal generator of an αresolvent family S α (t) which is compact for t ∈ (0, T ]. For T > 0, it follows from [1,Theorem 3.4] that the problem (6)-(7) admits a solution.

Sliding control with same diffusivity
In this section, the system (1)-(4) is with A(x) = a(x) I n×n .

Backstepping mapping
We use a backstepping mapping with a kernel matrix function K = [k ij ] ∈ R n×n , i, j = 1, .., n, to convert the underlying closed-loop system (1)-(4) into the target one where ∈ R n×n will be chosen to ensure the asymptotic stability of system (10)- (13).
With the same procedure of backstepping design as the one in [24], one can derive the kernel matrix function K(x, y). First we insert (12) and (3) into the firstorder spatial derivative of (9) (with x = 0) to yield since W (0, t) = U (0, t). We next calculate the Caputo time fractional derivative and second-order spatial derivative of (9) with (1)-(2), (10), (14) After this, we obtain that K(x, y) is given by The main critical feature of (15)- (19) is in the presence of relation (19). In order to simply to solve the above kernel PDE, we choose Ψ 1 = ψ 1 I n×n . Then we obtain that Ψ 1 K(x, y) ≡ K(x, y)Ψ 1 . With this, (15)- (19) reduces to (15)- (18). The following result is in order.
Proof. This proof is straightforward by using the one in [9, Appendix A].
Remark 2. For a more general Ψ 1 ∈ R n×n , the wellposedness of (15)- (19) remains an open problem, which could be interesting and will be considered for the future topic.
We next show the invertibility of the transformation (9), whose inverse transformation is given by with H = [H ij ] ∈ R n×n , i, j = 1, .., n. By the same implementation as those made for (15)-(18), we have Obviously, this kernel matrix PDE has a very similar structure to (15)- (18). In this sense we can obtain the existence of the inverse kernel by applying the same approach. The result is summarized as below.
Lemma 4. Suppose that a ′ (0) ≤ 0. The kernel matrix PDE (21)-(24) possesses a unique solution, and this solution is bounded and twice continuously differentiable in T .
Remark 3. We observer that K and H kernel equations have the similar feature to the one of same diffusivity in [9, Section III]. Herein they are extended to use for designing backstepping sliding mode controllers to cope with input disturbances.

Control strategy
In order to obtain the control input V c (t), we take the space derivative of the transformation (9) and insert this derivative, (9) into (3), (12). Following this issue, we have the control input We next employ the sliding mode algorithm to design U c (t) as follows where U 1 , U 2 ∈ R n×n , signW (1, t) = (signw 1 (1, t), ..., signw n (1, t)) T ∈ R n , U 1 and U 2 are subject to appropriate inequalities that will be derived later. In addition, sign · stands for the multi-valued function: Relation (26) depends on W (1, t) which is given by (10)- (13). And the discontinuous boundary controller (26) will enable to make the target system (10)-(13) asymptotically stable. Due to the invertibility of the transformation (9), the original system (1)-(4) will be also asymptotically stable under the controller (25).
Since the control input (26) contains discontinuous right hand side, W (x, t) is a solution of the target system (10)-(13) driven by this discontinuous controller, in the sense of Filippov [15]. In this case, the solution of the original system (1)-(4) with the discontinuous controller (25)- (27), is specified in the sense of Filippov. Please refer to [6,21] for more details on extension of the Filippov concept towards the fractional order systems.

Well posedness of closed-loop systems
To implement the proposed control strategy, the main question is if the original system (1)-(4) does indeed have a solution, or in other words, if the target system (10)-(13) does possess a solution (due to the invertibility of (9)). The next argument will answer this question. Consider the target system (10)- (13). Similar to the integer order PDE boundary control problem in [11,Theorem 3.3.3], we introduce the variable to transform (10)- (13) into an equivalent one with discontinuous right-hand side. For this discontinuous Wsystem, an operator A is introduced as follows: with the domain With A, this discontinuous system can be recorded as the fractional evolution equation in where ., e q s n (1−x) ), and ǫ(t) is the initial value based on the above variable.
In accordance with [3, Section 3.2], we split the state space H 1 into two subspaces S 1 and S 2 by a hypersurface Σ such that In this sense they are implicitly characterized by In (30), the right hand side P(t, η 1 (τ (t)), η(t)) is continuous or smooth on S 1 , S 2 respectively, however, it is discontinuous across Σ. A most important way to address this problem is to consider the multivalued extension F (t, η 1 (τ (t)), η(t)) as follows: where co{f 1 , Proof. Returning to the controller (26), for any T > 0, T ∈ (0, h 0 ], W (1, t) is continuously differentiable on [0, T ] due to the Caputo fractional derivative definition (Definition 1), which implies that the right hand side of (26) is locally bounded. This, together with (5), implies that F is locally bounded on t ∈ [0, h 0 ]. With a(·) ≥ a min > 0 and the constant matrix C, we further get that A is the infinitesimal generator of an analytic C 0semigroup. Finally, from Lemma 1 we obtain that the problem (32)-(33) has a solution on t ∈ [0, h 0 ].
Returning to the target system (10)-(13), we de- and the initial space for the target system Remark 4. Differently from the integer order nondelayed system [29], the resulting target system (10)-(13) is a coupled fractional PDE system and has a solution in the sense of Filippov, which can be viewed as the generalization of the Filippov method into the fractional order systems. And the solution of (10)- (13) is continuous. In [29], the well posedness of the integer order system has been mentioned (could be verified in accordance with [11, Theorem 3.3.3]) without rigorous demonstration.
Then by the invertibility of (9), we obtain the existence of the solution of the original system (1)-(4). Due to space limitation, the present work focuses on the existence of the solution whereas the rigorous demonstration of its uniqueness remains an open problem which is beyond the scope of this work. The uniqueness of (30)- (31) in question under the assumptions, is actually verified in accordance with [21,Theorem 6.14]. With this, we have the uniqueness of (10)-(13) (and, therefore, of (1)-(4) as well). In this case, the well-posedness of (1)-(4) is obtained.

Asymptotic stability analysis
In this section, we will derive the asymptotic stability of the delayed system of coupled fractional PDEs (1)-(4) via the fractional Halanay's inequality and the fractional Lyapunov method.
Proposition 2. Suppose that the target system parameters and controller parameters satisfy the following condition where q s min = min{q s 1 , .., q s n }, p s min = min{p s 1 , .., p s n }. Then for any J 1 ∈ Z 1 the target system (10)-(13) with the boundary sliding mode controller (26)-(27) is asymptotically stable.
Proof. The proof is straightforward since the original system (1)-(4) is related to the target system W by the invertible transformation (9) and the W -system is asymptotically stable (see Proposition 2).
4 Sliding mode control with distinct diffusivity

Deviation of kernel PDE and design of sliding mode control
In this section, we consider the distinct diffusivity case of A(x) = diag(a 1 (x), ..., a n (x)), a i (x) > 0, i = 1, ..., n. By the transformation (9) with a kernel matrix K = diag(k 1 , ..., k n ), we convert the system (1) ψ n−1,1 ψ n−1,2 · · · · · · ψ n−1,n ψ n1 ψ n2 · · · · · · ψ nn (x) into the following tar- In a similar fashion, we have Here we pick Ψ 1 = diag(ψ 1 1 , ..., ψ n 1 ) for simplification, we have Ψ 1 K(x, y) ≡ K(x, y)Ψ 1 . With this and A(x) = diag(a 1 (x), ..., a n (x)), (52)-(56) can be recorded as Remark 5. For the disturbance-free control problem of FRD system with distinct diffusivity in [10], the kernel equation is simplified by assuming a diagonal structure of the kernel matrix. There, the last term in the diagonal matrix is zero for designing underactuated controllers. In this work, the diagonal kernel matrix does not include the zero term and is used to design backstepping controller for the delayed system in presence of perturbations. In addition, if Ψ 1 and K(x, y) ∈ R n×n take the general expressions, then the analysis of wellposedness of (52)-(56) will need more techniques and resort to the one in [35].
The well posedness of PDE (57)-(60) has been proved (see [8, Lemma 1]). In order to prove our below results, the bound of k i (x, y) is needed. In this case, the mathematical computation procedure is shown as follows.

Using a transformation
we have the following integral equation of (57)-(60) By the similar procedure as the one presented in [10, Section 3.2], it follows by the assumption This implies that a i (0) , i = 1, 2, · · · , n.
Consider the inverse transformation (20) with H = diag (h 1 , ..., h n ). By a similar argument in [10, Section 3.2], we obtain that h i is given by for i ∈ [1, n]. Then we also obtain the bound of h i (x, y) as below The structure of the h i -kernel PDE (62)-(65) is very similar to the one of k i -kernel PDE. With the same approach, we can immediately have that the PDE (62)-(65) is also well posed. Therefore, the target system (47)-(50) can be converted into the system (1)-(4) by the transformation (20) with H = diag(h 1 , ..., h n ). For the distinct diffusivity case, we also choose the sliding mode control strategy (26)- (27) for U c (t) and the combined backstepping/sliding mode controller (25)- (27) for V c (t) to stabilize the target system and the original system respectively.

Well posedness of closed-loop system
In the distinct diffusivity case, since the transformation (9) is invertible, the well posedness of the original system (1)-(4) is equivalent to the one of the target system (47)-(50). In this sense the problem here is to show that the target system (47)-(50) is well-posed. Similar to the argument in Section 3.1, we use the variable (28) to transformation the W -system (47)-(50) into a discontinuous system. For this W -system, we also introduce the operator A with the domain D(A) in the Hilbert space H 1 . With the operator A at hand, we can rewrite this W -system as the fractional evolution equation in (30)- (31). For this distinct diffusivity case, the expression of the fractional evolution equation is same as the one of the same diffusivity case expect dξ dy. We can use the same analysis as the one of the same diffusivity case in Section 3.1 to obtain the following result for the system (47)-(50). Proposition 3. Consider the initial space Z 1 for the target system (47)-(50). There exists a Filippov solution of (47)-(50) for t ≥ 0 with the IC J 2 ∈ Z 1 .

Stability analysis
Before continuing, we provide a useful lemma which will be used for the main result.
Proposition 4. Suppose that the target system parameters and controller parameters are chosen as follows where G = max Consider the initial space Z 1 for the target system (47)-(50). For any J 2 ∈ Z 1 the system (47)-(50) with (26)- (27) is asymptotically stable.

Simulation study
For the illustration of the presented boundary sliding mode control approach, we consider a coupled FRD system with time varying delay which is inspired by com-bustion processes treated in [36]. In order to demonstrate the effective of the approach, we test the capabilities of the proposed boundary synthesis in simulation study. Due to space limitation, the boundary sliding control problem of two coupled fractional PDEs with the same diffusivity is treated, whose system is subject to input disturbances, acting on the boundary at x = 1 side. The behavior of the coupled system for the distinct diffusivity case is completely analogous.
To show the performance of the sliding mode algorithm, we first give the behavior of the plant (1)-(4) with only the backstepping controller which is depicted in Fig. 3. In comparison, the response of the closed-loop plant (1)-(4) with the whole controller (25)- (27) is displayed in Fig. 4, which shows the sliding control law U c (t) has suppressed the perturbation Φ(t). In Fig. 5, the state norm of this closed-loop plant is shown. This, together with Fig. 4, gives that the closed loop trajectory asymptotically converges to zero. Finally, Fig. 6 describes the average components of v 1av (t) and v 2av (t) of the discontinuous controllers (control input) v 1 (t) and v 2 (t), which is extracted through the low-pass filtering to estimate the reversed disturbance −ψ 1 (t) and −ψ 2 (t). While passing through the low-pass filtering, the discontinuous signal is smoothed out, and the chattering behavior unexpected in practice, is thus attenuated.

Concluding remark
In this work, we consider a time varying delay and a disturbance at the boundary control input for a coupled time FRD system where the coupling occurs in the system interior. By developing the combined backstepping/sliding mode control strategy, the perturbation has been rejected and the boundary feedback stabilization problem has been solved. An interesting research topic for future work is to take the state feedback stabilization and output feedback stabilization problems into account, by combing the backstepping method and the sliding mode control algorithm for fully coupled perturbed time fractional PDEs (with general kernel PDEs). Data availability My paper has no associated data.
Code availability Code will be made available on reasonable request after the publication of this paper. Authors' contributions Juan Chen conceived the original idea and prepared the initial draft, Bo Zhuang contributed to the code and software, YangQuan Chen conducted the methodology and revision, and Yajuan Yu supported validation.

Compliance with ethical standards
Conflict of interest The authors declare that there is no conflict of interest regarding the publication of this paper.  . 6 The average components of v 1av (t), v 2av (t) and the sign-versed disturbances −ψ 1 (t), −ψ 2 (t).