Leader-Follower Formation Tracking for Differential-Drive Wheeled Mobile Robots with Uncertainties and Disturbances based on Immune Fuzzy Quasi-Sliding Mode Control

Mobile robotics has drawn attention of researchers due to its variety of applications. In this paper, for nonholonomic diﬀerential-drive wheeled mobile robots with an embedded PD dynamic controller, a robust kinematic controller based on sliding mode theory is designed considering the trajectory tracking problem and its extension to the formation control by multiple robots using the leader-follower separation-bearing strategy. A classical sliding mode control has some drawbacks such as the chattering phenomenon and the requirement of a priori knowledge of the disturbances eﬀects boundary for its design. Thus, this paper aims to propose the design of an approach inspired by an immune regulation mechanism, also using a simple fuzzy boundary layer method, in replacement of the discontinuous portion of the controller, adjusting its gains with no need to know the uncertainties and disturbances boundaries as well as suppressing the chattering. The eﬀectiveness of the proposed controller is shown by simulation results and its stability is proved by the Lyapunov theorem.


Introduction
Wheeled mobile robots (WMRs) have applications in several areas, whether in industrial tasks or not, being attractive the possibility to design them to move without human intervention and most of them are already available on the market [1,2].Studies on multiple robots systems have grown due their better spatial distribution capacity, adaptability and scalability to perform varied tasks compared to a single robot [3][4][5].The differential-drive wheeled mobile robot (DWMR) is a typical example of a nonholonomic system and, despite the movement constraints, it is characterized by a simple and low-cost mechanical structure, widely applied in control experiments [6,7].A control problem in mobile robotics is the trajectory tracking, in which the robot must track a feasible trajectory with an associated temporal law.Whereas the formation control problem consists of driving each robot of a set to enter and maintain in a specific formation.The most popular techniques to ensure this are [8][9][10]: generalized coordinates, virtual structure, behavioralbased, leader-follower separation-separation and leader-follower separation-bearing approaches.In this paper, the last one is used for simplicity and the fact that its disadvantage associated to direct error propagation does not significantly affect formations with few DWMRs [3].Decentralized structure and homogeneous architecture (each robot has its controller and the same control hardware, respectively) are adopted.This strategy can be intuitively interpreted that a DWMR, as the leader, tracks some reference, while other DWMRs, as followers, track the leader with a separation distance and bearing angle [11,12].
This paper describes a simple and robust kinematic controller denominated first order sliding mode control (here simply denoted by SMC), depicted by a high-speed switching control law, whose task is to drive the trajectory of states onto a specified region of state-space (known as sliding surface) and keeping it there (property called by sliding mode) [13].It is often desirable to obtain means of developing control systems that make the minimum use of energy or that minimize the control effort (which may be related to actuator limits) without loss of performance, being an interesting concept to take inspiration from nature to design them [14,15].The artificial immune system (AIS) is an example of a computational intelligence bio-inspired approach.Typically, rather than replicating the biological immune system (BIS) as a whole, specific immune functions or models are considered in the design of the AIS [16].In mobile robotics, for trajectory tracking problem, most works propose algorithms based on the clonal selection theory for parameter optimization purpose, as in [17][18][19][20][21][22].Among these works only the last one considers the sliding mode control, but the formation control is not treated.When considering the SMC, most researches propose control structures inspired on the adaptive humoral/feedback immunity model, applying to different systems and problems, as in [23][24][25][26][27]. But, these works do not consider DWMRs.This type of AIS scheme is widely applied in PID-type controllers gain adjustment as in [28] for a WMR control (but without considering disturbances effects), and [29][30][31][32][33][34][35][36][37], for other control systems, to mention a few.In SMC, since the uncertainties are usually unknown, difficult to model accurately and rarely estimated a priori, considering a conservative control law with a large switching gain is often necessary to guarantee the system stability and a fast tracking convergence, but it can cause an undesired great amount of high-magnitude chattering, actuator saturation and unnecessary energy loss [13,[38][39][40].For this reason, an AIS inspired by the humoral immune response mechanism is proposed to perform an adaptive adjustment of the magnitude of the control effort by replacing the constant gain of the switching portion of a conventional SMC.To improve performance and chattering mitigation, a fuzzy boundary layer approach is also used.Differently from related works, for an online adjustment of the non-linear function that represents the reaction effect of immune cells in the AIS, a novel set of fuzzy rules substantially based on the principles of the humoral immune response mechanism is proposed for the leaderfollower formation kinematic controller applied to DWMRs under uncertainties and disturbances incidence, which are characterized by a non-linear, nonholonomic and under-actuated system, a relevant challenge in control field.Avoiding the use of time delays and parallel control is also concerned.So, the contributions of this work are: It is important to point out that this paper describes the problem of integration of the AIS-SMC controller (kinematic control) with a PD control (dynamic control) as well.Similarly to [41], velocity control signals are generated by the kinematic control to compensate posture tracking errors providing robustness and are also set as references for the PD control, which generates forces and torques that actuate on the DWMRs to compensate auxiliary tracking velocity errors.
This paper is organized as follows: Section 2 presents the kinematic and dynamic models of the DWMR, the trajectory tracking and formation control problems are described as well.In the Section 3, the kinematic and dynamic control structures (SMC and PD control) are presented.The Section 4 is dedicated to perform the control design considering chattering attenuation.The Section 5 describes the separation-bearing modeling.The simulation results are shown in Section 6. Conclusions, in Section 7, end this study.

Problem Formulation
The leader DWMR is indicated by the subscript i and the follower DWMR, by j.When the equations refer to the leader only or about DWMRs modeling in general, without loss of generality, subscript i will be omitted for readability.

Kinematic and Dynamic Models of the DWMR
A suitable mathematical model that represents the main kinematics and dynamics, namely posture model, is widely treated in literature [3,[41][42][43][44][45].A posture model in state-space of DWMRs with uncertainties is described by [3,43,46]: (1) with where q = [x y θ] T is the posture vector, with the DWMR Cartesian localization (x, y), and orientation T is the DWMR velocity vector; τ = [τ r τ l ] T is the torque vector, with torque in the right wheel, τ r , and left wheel, τ l ; m is the DWMR's mass; I is the moment of inertia; d is the distance between the center of mass and the intersection point of the symmetry axis of the robot and the wheel axis; b is the distance from the center of an actuated wheel to the DWMR symmetry axis; and r is the radius of the actuated wheels; M (q) and H(q, q) are known smooth nominal functions associated to the inertia and the centripetal and Coriolis matrices, respectively; and E(q) is an input transformation matrix.Representing uncertainties related to imperfect knowledge of physical parameters, actuator and sensor unmodeled dynamics, and external disturbances, the term τ b can be as [43]: The DWMR posture model, Eq. ( 1), presents the kinematics and dynamics in the first and second line respectively, pointing out that it will be used as the basis for the controllers synthesis only.

Trajectory Tracking
The reference trajectory is generated by a virtual DWMR with the reference posture and velocities known.The aim is to find a smooth control velocity v c such that lim t→∞ (q e ) = 0, where q e is the posture error and v r is the reference velocity vector.And then the control torque for Eq. ( 1) is computed such that v → v c as t → ∞ [47].

Leader-Follower Formation Control
The aim of separation-bearing formation control is to find a velocity control input v cj for the follower DWMR j such that [48]: where L ij and Ψ ij are the measured separation distance and bearing angle of the follower DWMR j to leader DWMR i, with L ij d and Ψ ij d being desired separation distance and bearing angle respectively.In addition to meeting these conditions, Eq. ( 4), this same smooth velocity input v cj must also satisfy lim t→∞ (q rj − q j ) = 0. Then torque τ j for dynamics, Eq. (1), is computed so that lim t→∞ (v cj − v j ) = 0. Achieving this every leader i and follower j ensures that the entire DWMRs set tracks the formation trajectory [3].

Dynamic Control
The dynamic control aim is to compensate known torques and forces, to ensure fast-tracking of the auxiliary velocity tracking errors, v e .For the PD control design, as [41,43,49], one considers: applied to the system, Eq. (1).To achieve fast convergence of v c , the vector u = [u v u ω ] T is a control input that will be designed as the PD control, being generated in frequency domain by , with the derivative filter parameter gains N v and N ω being positive and adjusted to achieve stability with good time response performance in spite of the neglected dynamics, avoiding chattering [43].
The stability for the PD control law, Eq. ( 5), is verified by the candidate function [3,49,50]: which is positive definite.With v e = v c −v and the proportional gains matrix k p = [k pv 0; 0 k pω ], and the derivative gains matrix k d = [k dv 0; 0 k dω ], considering the error dynamics as [3]: M (q) vc + H(q, q)v c = γ(q, q, v c , vc )ϕ ; M (q) ve = γ(q, q, v c , vc )ϕ − H(q, q)v e − E(q)τ , the time-derivative of Eq. ( 6) is obtained as: λ min {k p } is the minimum singular value of k p .VP D < 0 while the term in parentheses is positive.This dynamic control architecture ensures fast convergence of auxiliary velocities tracking errors, as seen in [41,43,49] and it is proved asymptotically stable by Lyapunov's theory considering a constant reference velocity.However, in this paper, the reference velocities are time-varying, which implies the stability of the closed-loop system loses its asymptotic feature.For this reason, the residual error will be handled by the SMC at the kinematic control.It is also noteworthy that the kinematic control could handle the entire control problem without the PD controller, but since the considered PowerBot DWMR is a closed platform with an internal dynamic control architecture, the use of PD control is maintained.

Kinematic Control -SMC
The aim of leader kinematic control is to determine the velocities required so that the leader DWMR i can track a feasible trajectory which can be described by a virtual reference DWMR.Being q r = [x r y r θ r ] T and v r = [v lr ω ar ] T as the reference posture and velocity vectors respectively, the virtual DWMR kinematics is modeled as: Based on the reference and DWMR postures, the posture tracking errors are obtained in the inertial frame to the DWMR frame, and then the leader DWMR i posture error is derived as: The error dynamics of the closed-loop system for trajectory tracking, using Eqs.( 1) and ( 2), after algebraic manipulations, is given by: qe Due to nonholonomic constraint [11] and d ̸ = 0, the leader DWMR i orientation will not be equal to the virtual DWMR orientation while making a turn.For this reason, the reference orientation cannot be chosen as θ r = θ.
Suppose that it is chosen: The matrix J * r (q * r ) is invertible for d ̸ = 0. Thus, seeing Eq. ( 9), auxiliary reference velocities (representing the desired movement in the task space) can be obtained from the inverse kinematics by Taking into account the auxiliary reference angular velocity ω * ar from Eq. ( 13), based on a generic modeling of nonlinear systems adding the vector d b to consider disturbances incidence, the error dynamics, Eq. ( 11), can be rewritten as: where A 0 (q e , t) and B 0 (q e , t) refers to the vector and matrix of nominal parameters, respectively, and v c (t) is the vector of control inputs.Since the "perfect velocity tracking" (v = v c ) does not hold into practice, some bounded auxiliary velocity tracking errors v e are generated, which can be transcribed as uncertainties and disturbances effects d0 for the kinematics [41,47].The disturbances effects d0v and d0ω are assumed to be bounded by positive constants as For SMC control synthesis, the first step is the choice of the sliding surfaces, which are, in general, defined by linear surfaces of the control errors and its vector must have the same dimension as control velocities vector v c , so σ v and σ ω are related to v lc and ω ac , respectively [43].Thus, based on Eq. ( 14), the sliding surfaces are selected as [3]: where Λ 1 , Λ 2 and Λ 3 are positive constants.
However, to solve a nonsingularity problem of B 0σ , as described in more detail in [3,41,45], new sliding surfaces are required to be chosen as: It is pointed out that the aim of the SMC is to obtain a sliding motion restricted to σ(q e , t) → 0, σ(q e , t) → 0, σ * (q e , t) → 0 and σ * (q e , t) → 0.
For the stability analysis, the Lyapunov candidate function is chosen as: which is positive definite (V > 0).To ensure the attractiveness of the sliding surfaces, the control law is derived as: where G and P are positive definite diagonal matrices, sign(σ * is a discontinuous function given by sign(σ Seeing Eq. ( 16) and the control law Eq.( 19), the derivative of Eq. ( 18) results in: Since σ * T P σ * ≥ 0, the condition V ≤ 0 is satisfied when λ min {G} > d0 M , being λ min {G} the minimum singular value of G and d0 M , the disturbances maximum effect.Hence, it is possible to affirm that if the uncertainties and disturbances are better estimated, the results will be better.

Leader Control Design
Considering Chattering Attenuation

Humoral Immune Response Mechanism
The main role of the BIS is to identify and eliminate foreign bodies called antigens, being the collective and coordinated response of the immune cells known as immune response.There are two types of immunity, which work together to carry out the organism's defense [51][52][53]: the innate immune system, which is the first line of defense, reacting to all pathogens by immediately destroying them; and the adaptive immune system, which is directed against specific antigens, and whose immune responses are classified into cell-mediated immunity and humoral immunity.This paper considers the last one to design the proposed AIS control structure.The cells responsible for adaptive humoral immunity are white blood cells (leucocytes) called lymphocytes, being the main ones B cells and T cells.The B cells are lymphocytes associated with the production of antibodies (Ab) that attack a specific antigen (Ag).Whereas T cells, to be activated, need antigen-presenting cells (AP C).The immunity regulation is mainly performed by T helper cells (T h ), which secrete interleukin IL + that stimulate B cells to produce antibodies, and suppressor T cells (T s ), which, by producing interleukin IL − , restrict the action of T h and B cells in order to prevent excessive immune responses, leading to the response stabilization.The simplified humoral immunity behavior can be illustrated according to the Fig. 1.Generally, projects of AIS based on this type mechanism focus on the restrictive action on B cells [31].So, similarly to [30], the corresponding parts of humoral-type AIS structure can be obtained by the following equations: where η = k s (ε)/k h (ε) represents a proportionality factor describing the interaction between T h and T s cells.The meaning of each AIS term is shown in Table 1.And the Table 2 outlines the control system course in comparison to immune regulation process [37,54,55].The terms k h and k s can be variables based on the antigens amount (e.g., a linear/nonlinear function or an adaptive parameter), as in [24], [25] and [27].In this work, they are considered as positive constants, as in [23] and [26].

Immune Fuzzy Quasi Sliding
Mode Control (IFQSMC) Given the concepts described, the immune systems inspired approach is adopted to minimize the control effort and compensate for uncertainties and Nonlinear function that represents the immune reaction effect.
Table 2 Immune reaction and adaptive control.

Immune reaction Adaptive control
Facing an invasion of foreign materials, when the first line of defense (innate immunity) fails, the amount of antigens is higher and the production of antibodies is expected to increase rapidly, so the inhibitory action of Ts cells must be contained to promote greater stimulation.
If robustness fails (control action is unsatisfactory or system is unstable or disturbed), an estimate of the limit of disturbances is presented to generate an "antibody" equivalent to the modification of the control system.With antigen reduction, the amount of antibodies is not expected to increase continuously.In other words, the containment on Ts cells should decrease.
Based on some performance measure, the parameters of the control system can be adaptively adjusted.
When the most of the antigens have been removed, the suppression stimulus of Ts cells is expected to increase rapidly to restrict B cells and antibody production, leading to stabilization of the immune reactions.

When
the observed errors present an indicative of reduction, the control system output is computed in a such way that minimize control efforts and achieve control system stability.disturbances adaptively.For this, we established that the antigens amount is the sliding surface σ * , as in [23][24][25][26][27], and total stimulation received by B cells for antibody generation corresponds to the sliding mode portion gain G AIS of the AIS-SMC controller.In this way, by Eqs. ( 23), ( 24) and (25), the nth AIS output signal is denoted as: with k hn and k sn being positive constants.
The proposed AIS control structure for the gain G AIS is shown in Fig. 2.This immune regulation inspired mechanism can smooth control effort magnitude, but the discontinuity still remains due to the signum function.Thus, to improve chattering attenuation, it is introduced a boundary layer in the control structure.For it, sign(σ * ) is replaced by the saturation function, sat(σ * ), resulting in the IFQSMC.
where the nth element of sat(σ * ) is defined as: in which β n > 0 defines the nth value of the boundary layer thickness.This alternative has a price that the attractiveness is no longer ensured to zero but to a closed region defined by the boundary layer [56].Despite this, it is able to smooth the discontinuous control signal while the states trajectory is in zero vicinity.
It can also be seen, from concepts summarized in Table 2, that the nonlinear function f (•) plays an important role related to regulation of the stimulus received by B cells, since the reaction can be defined by the activation or containment of the effect of T s cells.In some works, f (•) is described by an exponential function based on antigen consistency as in [24,25,27,31,32,34].A fuzzy inference system is also used as in [23,26,29,30,[32][33][34][35][36][37], which usually consider the total stimulation received by B cells and its variation with a time delay of immune response.In this paper, in order to also provide greater adaptability and without the need to apply a time delay, we investigated the reaction effect based on the antigen amount and its variation, according to the principles described in Section 4.1.
So, taking advantage of the fuzzy systems universal approximation property, simplicity and intuitive linguistic interpretation, a Mamdani fuzzy inference system (Mamdani FIS) is designed to compute the nonlinear function f (σ * , σ * ).For the FIS, with N input variables and M output variables, triangular membership functions was chosen, and it was adopted the singleton fuzzification, the product inference engine, the center average defuzzification method and a fuzzy rules base composed of K "if...then..." rules, in which the k th rule can be denoted as m are labels of fuzzy sets in input and output spaces, respectively.Adopting this structure, the FIS output is obtained by [57]: where c (k) m is the center value of the mth output membership function (MF) in the k th fuzzy rule, ) is the product of MFs' mappings for all input variables in the k th fuzzy rule, being µ n the membership function of the input x n , the signal y m (m = 1, . . ., M ) is the mth crisp output of the FIS, and ξ(x) = ξ (1) In this paper, the non-linear function f (•) is the output of the fuzzy controller, whereas the sliding surface σ * and its derivative σ * are the inputs.The fuzzy rules were designed according to the immune regulation mechanism principles summarized in the Table 2 and can be inferred as shown in Table 3, where NB is Negative Big, NM is Negative Middle, NS is Negative Small, ZO is Zero, PS is Positive Small, PM is Positive Middle and PB is Positive Big.Since σ * can assume positive and negative values, it can be interpreted, in summary, that there is an antigen amount in distinct semi-cycles: positive and negative.This means that a positive value for σ * and a positive variation ( σ * > 0), or a negative value for σ * and σ * < 0 denote an increase in antigen amount, whereas σ * > 0 and σ * < 0, or σ * < 0 and σ * > 0 indicates a reduction of antigens.If σ * and σ * are close to zero respectively indicate that there is no significant amount and no significant variation has occurred.Moreover, f ≤ 0 represents a containment of the suppression effect, that is, the stimulus is enhanced, whereas f > 0 denotes a promotion of the inhibitory effect to smooth the control effort.Tuning the boundary layer thickness may also be convenient.In [58], a chattering measure is introduced as the absolute value of the derivative of the control signal, | vc |, considered as a FIS input variable with the sliding surface to compute the boundary layer alteration.With a exponential relation to regulate the control effort of the reaching law, this concept also inspired the conception of the immune based SMC in [27].Although this approach suggests a more dynamic adaptability by taking into account the chattering level, and since the immune control already relies on a regulation process from the f function by a FIS, to provide an online adjustment of the BL thickness, also avoiding time delays, we opted to implement a Mamdani FIS considering as input only the absolute value of sliding surfaces |σ * n | with the boundary layer thickness β n as the output variables, similarly to [59].
For the FIS design, the following principles can be brought up [27,59]: When the sliding surface is far from the sliding manifold σ * n = 0, a small boundary layer is required for a faster convergence; And, when it is closer to the sliding manifold σ * n = 0, a larger boundary layer is desired to mitigate chattering.According to these concepts, the fuzzy rules is designed as presented in Table 4, where VB means Very Big, B is Big, M is Middle, S is Small and VS is Very Small.

Rule
Antecedent Consequent The fuzzy output is obtained similarly to the Eq. ( 29).Just note that while, in the FIS designed for AIS function f (•), each output is related to two inputs (i.e., f v to σ * v and σ * v as well as f ω is to σ * ω and σ * ω ), this fuzzy system used to compute the boundary layer values have each output associated to only input (i.e., β v to |σ * v | and β ω to |σ * ω |).Thus, the sliding mode portion of the IFQSMC with the fuzzy boundary layer approach can be established as depicted in Fig. 3
The stability analysis is performed choosing the Lyapunov candidate function similarly to Eq. (18).When considering the case where σ * > β, since sat(σ * ) = sign(σ * ), one has: The system trajectory convergence to the boundary layer is guaranteed by the conditions established by Eq. (30), which is similar to Eq. ( 22), so the same conclusions on the stability analysis are valid, being the boundedness of each diagonal element of G AIS , i.e., g n AIS , is ensured by considering that the gains k hn and k sn are positive constant values as well as that function f n is bounded, since it is estimated by a Mamdani FIS, and can assume negative and positive values such that −f − n ≤ f n ≤ f + n .Note that if k sn = 0 or f n = 0, the IFQSMC is the same as a traditional SMC where g n AIS = k hn .Since g n AIS > 0, then k hn − k sn f n > 0, and thus, having established the gains, f n need to be upper bounded by f n < k hn /k sn .Therefore, it is important to take these conditions into account in the FIS design (for this reason, in this case, f was delimited by a "positive middle" value, in such a way that it meets these analysis thoughts).
Note also that k sn f n < 0 represents a restriction on the suppression effect, promoting the increase of the gain g n AIS , while 0 < k sn f n < k hn indicates a activation of the inhibitory effect, providing a decrease of the gain g n AIS in order to minimize the control effort.
The control within the boundary layer changes to a continuous and high-gain signal (sat(σ * ) = σ * /β).Thus, in this case, selecting the control law, Eq. ( 27), V results in: Thus, V < 0 is satisfied when: suggesting that there is a lower limit value of σ * n at which the inequality is satisfied.Note that a wider boundary layer can influence an error increase and higher gains will have to be.Thus, it can be estimated that the errors decreasing is related to suitable values of gain g n AIS and boundary layer β n .Moreover, since it is computed by a Mamdani FIS, β n is bounded such as 0 ≤ β nm ≤ β n ≤ β n M .

Leader-Follower Trajectory Tracking Control
The extension of the trajectory tracking problem to the formation control is done through the leader-follower separation-bearing strategy, whose aim is to find control velocities in such a way that it meets the conditions in Eq. ( 4), i.e., that the separation distance errors (L e = L ijd − L ij ) and the bearing angle errors (Ψ e = Ψ ijd − Ψ ij ) tend to zero.It means that the aim of follower DWMR j kinematic control is to establish the linear and angular velocities needed to track the leader DWMR i with a given separation distance and a bearing angle.To avoid collisions between them, the distance L ij is measured from the center of wheels axis of the leader DWMR i to the center of the front of the follower DWMR j.So, the value of L ij and Ψ ij can be as: where The kinematic model of the follower DWMR j can be described as: where the follower DWMR j posture is given by: Since the follower DWMR j tracks the leader DWMR i, its reference posture vector can be obtained as follows: So, the posture error is described by: (38) Replacing Eqs. ( 36) and ( 37) into Eq.( 38), and considering θ ij = θ i − θ j , using trigonometric transformations and simplifications, one has: To obtain the error dynamics, it is necessary to calculate the derivatives of L ij and Ψ ij , being the desired values L ijd and Ψ ijd known constants.So, calculating the time derivatives of L ijx and L ijy : Considering Eq. ( 33) and Due to the separation-bearing formation control aim and to the non-slipping constraint [11], the orientations of each DWMR in the formation will not be equal while making a turn.For this reason the orientation of each DWMR cannot be chosen in such a way that θ jr = θ i .Thus, from the inverse kinematics, being θ * ej = θ * jr − θ j , the auxiliary reference angular velocity is obtained as: Describing qej by using generic modeling of nonlinear systems and considering the effect of uncertainties and disturbances, one has: From the posture error dynamics, Eq. ( 42), the sliding surfaces for the follower DWMR j is defined as: The control laws of the SMC and IFQSMC for the follower DWMR j are established similarly to Eqs. ( 19) and ( 27) respectively, emphasizing that the sliding surfaces are determined by Eqs. ( 43) and ( 44), as well as the stability analyses of the control systems are similar to the stability analyzes performed for the leader DWMR i. Considering h followers and proper assumptions as made in [3,11,48], the formation stability can be demonstrated utilizing the individual Lyapunov functions and their derivatives with the relative conditions being required and met, considering the closed-loop control systems for the leader DWMR i and the follower DWMR j [3].Thus, the following Lyapunov candidate function is chosen: where V i P D and V j P D are determined as Eq. ( 6) while V i and V j are established as Eq.(18).Since V j P D ≥ 0 and V j ≥ 0 for all j, and V i P D ≥ 0 and V i ≥ 0, thus, it can be conclude that V ij ≥ 0. Differentiating the Eq. ( 45), one obtains: where Vi P D and Vj P D are determined as Eq. ( 8) while Vi and Vj are established as Eq. ( 22).Since Vj P D ≤ 0 and Vj ≤ 0 for all j, and Vi P D ≤ 0 and Vi ≤ 0, thus, it can be conclude that Vij ≤ 0.

Simulation results
Simulation results using Matlab/Simulink software (version R2020b), with the total simulation time of 200s and the Dormand-Prince solver, are obtained for the trajectory tracking formation control, for the leader DWMR i and the follower DWMRs j, the PD control as dynamic controller, with control law as Eq. ( 5) (which can be described as τ = E(q) −1 (k d ve + k p v e ) in time domain [3]), integrated with the kinematic controllers SMC and IFQSMC, with control law established as Eqs.( 19) and ( 27) respectively.In the simulation scenario, it is considered a configuration model with the limitations of the PowerBot DWMR, which is a closed platform with internal PD controllers that tracks inputs variables v e and ω e with a sampling time of 5 ms.Furthermore, some modifications to the DWMR configuration model are necessary because the internal PD control signals generate generalized forces acting at the DWMR's center of inertia and must be converted into torques on the wheels to be taken as control inputs in the present configuration model.So, the signals of the PD controllers are premultiplied by the matrix T −1 v , which relates torques to the wheels with generalized torques, given as [43]: r , where 2b is width of the PowerBot DWMR.The specifications and parameters in nominal values of the PowerBot DWMR can be seen in [43] or in the manufacturer's manual.
The gains of the PD dynamic controller are adjusted as k pv = 40, k pω = 40, k dv = 20 and k dω = 20 as indicated by the manufacturer.The parameters N v and N ω are set to 10 to ensure good response time of the internal loop against the neglected dynamics and avoid chattering [43,49].
An eight-shape trajectory was adopted as reference trajectory and is formulated as [3]: The reference trajectory established by Eq. (47), in view of the acceleration and deceleration along it, has a complexity defined by the rapid changes in velocities, with linear velocity varying between 0.33 m/s and 1.05 m/s, and the angular velocity, between -0.18 rad/s and 0.18 rad/s.
To define the formation shape, the desired separation-bearing values L ijd and Ψ ijd are defined as shown in Table 5.And the Table 6 presents the initial conditions and the values of d of the DWMRs.The gains of the kinematic controllers are defined in Table 7, remembering that the positive constant gains g v and g ω are used only for SMC.The gains k h and η for IFQSMC are established in Table 8.For the fuzzy inference system used to compute the function f (•), the membership functions can be graphically represented by Figs. 4, 5, 6 and  7, being the parameters given in Tables 9 and 10, being the parameters given in Tables 9 and 10.
Similarly for the boundary layer thickness FIS, the MFs can be illustrated by Fig. 8 with the parameters given in Table 11, being the fuzzy discourse domain and fuzzy subsets of the inputs and the outputs assumed to be the same.Note that the same values of g v and g ω were adopted for k hv and k hω in order to verify the compensatory behavior of the IFQSMC.Moreover, an adequate choice of the parameters of the membership functions for the FIS output, mainly   referring to the its center values as established in the Eq. ( 29), is related to the chosen values for the terms of k h and k s , in such a way that the condition G AIS > 0 is obtained, i.e., f n < k hn /k sn , considering the Eq. ( 26) and according to the performed stability analysis.It is also suitable to emphasize that when considering the control for kinematics only, it is common to assume that the control velocities are equal to the DWMR velocities.However, this may not remain valid from a practical application point of view [47].By also considering the dynamics, the perfect velocity tracking assumption (v c = v) is no longer maintained, resulting in the generation of auxiliary velocity tracking errors (v e ), which can be transcribed as uncertainties and disturbances effects by the kinematics.So, in order to verify the compensatory characteristics behavior of the kinematic controllers in a perturbed case, simulations are performed considering load variations of mass and moment of inertia (in the time range of 40 to 160 s), as well as unmodeled dynamics introduced in the form of friction (during the entire simulation time).Additionally, external disturbances with unknown upper bounds are applied on the system at different time intervals (time-varying sinusoidal form at 3 to 20 s and 100 to 120s; constant signal at 50 to 70 s and 150 to 170 s).These uncertainties and disturbances are similar to those considered in [41].
The root mean square errors (RMS) is an evaluation method commonly used in the literature [41,43], so it is used as one of the ways to verify the controllers performances and is given by RM S(x) = 1 T x 2 1/2 , where T is the number of samples and x(k) is the value of the k -th sample.9,17).This leads to a larger RMS errors (see Table 12).Thus, as expected, this shows that, the effectiveness in terms of robustness against uncertainties and disturbances depends not only on an adequate choice of sliding surfaces but also of the gain values, and therefore a priori knowledge of the disturbances boundaries is fundamental to the sliding mode control design.Besides that the chattering is quite evident (see Fig. 18).

Simulation results for IFQSMC
By the results obtained in Figs.19 to 32, it can be seen that the IFQSMC was able to suppress the chattering as it is noticeable by comparing Figs.18 and 28.Satisfactory performance is obtained as can be seen by the RMS errors shown in Table 12.Moreover, the σ * -trajectory, i.e., the variation  of σ * with time, denotes a compact descriptor of the closed-loop behavior of the control activity, and can be a means of translating uncertainties and disturbances effects into observed errors.Similarly, the boundary layer thickness can describe the evolution of the control system.For this reason, it is particularly informative to plot σ * , +β and −β on a single diagram [60].So, for a better visualization of the behavior of σ * and β, graphics with zoom are presented as shown in Figs. 29 and  30, through which it is also noted that, in intervals without significant disturbances, σ * tends to remain within the boundary layer which enlarges to avoid chattering, whereas, in a perturbed case, its thickness is reduced in an attempt to improve

Conclusions
The aim of this paper was to propose a leader-follower separation-bearing formation control applied for DWMRs under incidence of uncertainties and disturbances in solving the trajectory tracking and formation control problems.
It is considered the traditional SMC, but this control need a proper choice of its gains to ensure the robustness and accuracy.The proposed kinematic controller IQFSMC, inspired on immune regulation mechanism, provided an online adaptability for the SMC gains, being able to not only smooth the chattering, but compensate uncertainties and disturbances transcribed as auxiliary velocity tracking errors derived from dynamics, without a prior knowledge of their upper bounds and avoiding unnecessary control effort as well.
Research can be extended to carry out experimental tests; generalization for other possible classes of mobile robots (i.e.unmanned aerial vehicles); and also the treatment of obstacle avoidance.
Intelligent motion control for four-wheeled holonomic mobile robots using fpga-based artificial immune system algorithm.
10 and 11), which ends up deteriorating the trajectory tracking and the formation (see Fig.
the σ * reach to it and enhance the tracking accuracy.The time-varying gains and parameters from the AIS also show that they are able to adapt to obtain the compensatory effect without exceed control efforts (Figs.31,32, 24, 25,16).Consequently, small errors are obtained, as well as the SMC and formation control aims are achieved(Figs.20, 21, 22, 23,27).

Table 1
Meaning of each AIS term.
h Action of T h cells.Ts Influence of Ts cells.S B Total stimulation received by B cells.

Table 3
Fuzzy rules for nonlinear function f (•)

Table 4
Fuzzy rules (adaptive variable boundary layer)

Table 5
Values of Ψ ijd and L ijd for followers.

Table 6
Initial posture and parameter d of the DWMRs.

Table 8
Gains k h and η for IFQSMC.

Table 9
Parameters of MFs for input σ * and σ * .

Table 10
Parameters of MFs for output f (•)

Table 11
Parameters of MFs for Fuzzy BL.

Table 12
RMS of the errors in a perturbed case.