Fractional-Order Sensing and Control: Embedding the Nonlinear Dynamics of Robot Manipulators into the Multidimensional Scaling Method

This paper studies the use of multidimensional scaling (MDS) to assess the performance of fractional-order variable structure controllers (VSCs). The test bed consisted of a revolute planar robotic manipulator. The fractional derivatives required by the VSC can be obtained either by adopting numerical real-time signal processing or by using adequate sensors exhibiting fractional dynamics. Integer (fractional) VCS and fractional (integer) sliding mode combinations with different design parameters were tested. Two performance indices based in the time and frequency domains were adopted to compare the system states. The MDS generated the loci of objects corresponding to the tested cases, and the patterns were interpreted as signatures of the system behavior. Numerical experiments illustrated the feasibility and effectiveness of the approach for assessing and visualizing VSC systems.


Introduction
Variable structure systems (VSSs) [1,2] have good feasibility and robustness. As concerns system theory, variable structure controllers (VSCs) are a relevant strategy when facing systems with complex dynamics [3][4][5]. VSCs apply switching control laws to alter the system dynamics and compel the system states to slide along a cross-section named the sliding surface. The system trajectories have two distinct periods, usually called the (i) reaching and (ii) sliding mode phases. In the reaching phase, the system states are forced toward a prespecified sliding surface in finite time. Once the system states reach the sliding surface, the sliding phase initiates, and the closed-loop system states slide toward the origin along the sliding surface. During the reaching phase, the invariance of the VSC is not guaranteed and the system response is quite sensitive to perturbations. In the sliding phase, the system response remains invariant for both parametric and nonparametric uncertainties [1,6].
Mechanical manipulators exhibit strong nonlinear dynamical effects that require appropriate control algorithms. When applying VSCs with mechanical manipulators [7][8][9][10], a common practice is to approximate the dynamics of each rigid link as a first-order linear model. The imperfect implementation of high-frequency switching in VSCs results in chattering at the control responses. This effect leads to a large stress on the actuation hardware and can excite vibrations in the structure. To reduce the two effects, several schemes have been proposed, involving some kind of smoothing, namely by changing the "on-off" switching algorithm, or by complementing the VSC with some adaptive or feedforward control actions [11][12][13]. However, the first-order reference model is not well suited to the system's intrinsic dynamics; therefore, a second-order model was proposed [2,14,15].
Fractional calculus (FC) extends the scope of the classical calculus to noninteger orders [16][17][18][19][20]. Fractional derivatives and integrals are helpful in control since they allow adapting existing integer algorithms to a fractional version with more degrees of freedom. Indeed, due its properties, FC emerged as a key tool in dynamics and control systems [21][22][23][24][25].
Embedding the FC concepts in VSCs has also been proposed [26][27][28][29][30]. Since we have the freedom of selecting the fractional order, we can take advantage of the new parameter to tune either the reference model or the control algorithm, affecting, thereby, the VSC switching action. Previous work revealed that combining FC and VSC leads to superior performance [31][32][33]. However, the algorithms require either numerical real-time signal differentiation or adequate sensors to obtain the fractional derivatives. We can cite recent works, namely Delavari et al. [28], who proposed a fractional-order sliding-mode controller to control a flexible-link manipulator, while determining the design parameters trough particle swarm optimization (PSO). Delavari and Heydarinejad [29] designed a fractionalorder backstepping-sliding-mode control for a class of fractional nonlinear systems with mismatched disturbances, which were estimated using a fractional nonlinear observer. Simulation examples showed the effectiveness of the control strategy. Wang et al. [34] investigated cable-driven manipulators under lumped uncertainties and proposed an adaptive fractional control scheme based on time-delay estimation. The controller included a timedelay estimation to compensate the unknown system dynamics, a fractional nonsingular terminal sliding-mode surface to ensure high precision in the steady phase, and a reaching law with an adaptive technique to obtain fast convergence, high precision, and reduced chatter. Simulations and experiments were performed to show the effectiveness of the scheme. Zhou et al. [35] proposed a deep-convolutional-neural network-based fractional terminal sliding-mode controller for rigid manipulators. The neural network compensated the uncertainties of the system. The chattering was mitigated, and the control strategy exhibited robust performance against uncertainties and disturbances. Ma et al. [31] proposed a quaternary fractional-order sliding-mode controller with fuzzy logic system, a neural network, and an adaptive law to control the teleoperated cyberphysical system. External disturbances, modeling uncertainties, and actuator faults were considered. Xie et al. [32] addressed a coupled fractional sliding-mode control, together with an obstacle-avoidance scheme, to control a four-wheeled steerable mobile robot. A modified near-time-optimal potential function was introduced to improve collision problems. Fuzzy rules and proper adaption gains were designed to mitigate chattering. Asymptotic stability and convergence were guaranteed for the closed-loop system. Delavari and Jokar [33] presented a fractional-order active fault-tolerant controller based on an adaptive nonlinear observer to detect, estimate, and compensate faults of a knee joint orthosis. The controller was based on fractional-order sliding-mode control, while the switching term was designed using fractional-order interval type-2 fuzzy logic. The strategy was proven to reduce modeling issues and chattering. Other examples can be found in [36][37][38].
Several researchers proposed a variety of fractional controllers both in the discretetime [39][40][41] and frequency [42][43][44] domains. The implementation of such controllers requires the use of approximation algorithms to obtain the fractional derivatives or the implementation of sensors exhibiting fractional dynamics. The construction of these sensors can take advantage of the advanced fabrication techniques of microelectromechanical systems (MEMSs) [45], allowing, therefore, directly obtaining the fractional derivatives [46].
We have well-known indices to assess the performance of controlled systems. However, in general, a single criterion is not sufficient to capture the dynamic details. Therefore, characterizing the system behavior is a multidimensional problem that can be tackled with recent computational tools [47,48]. Dimensionality reduction [49] plays a key role, since the data often exhibit a multidimensional nature. Dimensionality-reduction-based schemes try to preserve in lower-dimensional representations the information present in the original datasets. They include linear methods, such as classic multidimensional scaling (MDS) [50], principal component [51], canonical correlation [52], linear discriminant [53], and factor analysis [54], as well as nonlinear approaches, such as nonclassic MDS or Sammon's projection [55], isomaps [56], Laplacian eigenmaps [57], diffusion maps [58], t-distributed stochastic neighbor embedding [59], and uniform manifold approximation and projection [60]. Besides dimensionality reduction, these algorithms often allow direct information visualization.
Hereafter, we considered the MDS technique to evaluate and visualize the performance of fractional-order VSC. The test bed adopted consisted of a revolute planar robotic manipulator. The fractional derivatives required by the controller can be implemented either by numerical methods or by using sensors exhibiting fractional dynamics. Integer (fractional) VCS and fractional (integer) sliding mode combinations, with different design parameters, were studied. The system states for a number of test cases were obtained and compared by means of two alternative distance metrics, namely in time or in frequency. The information was input into the MDS. The algorithm then generated the loci of objects where each point corresponded to one test case. The objects patterns were interpreted as signatures of the system behavior. Several numerical experiments illustrated the effectiveness of the approach.
The paper has four sections. Section 2 addresses the mathematical background and key concepts useful in the remainder of the paper. Section 3 analyzes the dynamics of the VSC system by means of MDS, while adopting the integer VSC and fractional sliding mode (IVSC-FSM) and the fractional VSC and integer sliding mode (FVSC-ISM). The analysis is performed either in the time or frequency domain. Finally, Section 4 presents the main conclusions.

Fractional Integrodifferential Operators
The Grünwald-Letnikov (GL) fractional operator of order α ∈ R on y(t), denoted as GL a D α t y(t), is given by [61]: where [·] denotes the integer part, h is the time increment, and {t, a} ∈ R (t > a) are the upper and lower limits of the "differintegral" operation, respectively. In signal processing and control, the GL definition (1) can be approximated numerically using [62,63]: where T is the sampling period, L corresponds to the "memory length", and k are computed by [62]: The "memory length" L is often chosen considering: where Γ is the gamma function and δ 0 represents the maximum allowable error: In many practical applications, we consider a = 0 and adopt D α t to denote the generalized "differintegral" operator.
Alternatively, we can use the approximation: where Z {·} denotes the Z transform. Equation (6) gives the s → z Euler approximation, but we find often other algorithms in controller design, such as the so-called Tustin or trapezoidal scheme. The Euler and Tustin conversion techniques can be generalized in the scope of FC to: where Ψ 0 z −1 α and Ψ 1 z −1 α are approximation functions of zero and one order, respectively.
To obtain rational expressions, we need to truncate the Taylor series or Padé fraction. The two approximations Ψ 0 z −1 α and Ψ 1 z −1 α can be averaged with weights p and 1 − p, respectively, resulting in: For instance, the case p = 3 4 corresponds to the Al-Alaoui operator [64,65].

Variable Structure Control
When adopting the VSC in a manipulator, the k-th link (k = 1, · · · , K) is induced to mimic a first-order reference: with θ dk ,θ dk and θ k ,θ k denoting the desired and the true position and velocity of the k-th joint, respectively, σ k standing for the switching variable, and e k representing the position error. The expression s + λ k = 0 has an eigenvalue, λ k ∈ R, that characterizes the sliding phase dynamics. The VSC produces a control action T k (t) that forces the robot to mimic the reference (10). Often, the VSC follows an algorithm of the type: where the function sgn(·) returns the sign of its argument. When the VSC satisfies the condition: convergence is guaranteed.
In [66], it was observed that such a strategy leads to stringent requirements because the first-order model (10) allows discontinuous trajectories in the phase plane. However, mechanical manipulators have inertias that lead to the continuous evolution of the links' positions and velocities. Therefore, a first-order model requires high joint torques during transients. In practice, torques cannot be infinite, and thus, the phase plane trajectories are continuous. However, the demanding torque requirements saturate the robot actuators, resulting in a longer reaching phase, which is sensitive to perturbations. Therefore, a second-order model was proposed to mitigate these problems [14]: where ζ k denotes the coefficient of damping and ω nk stands for the undamped natural frequency. If the corresponding Laplace equation s 2 + 2ζ k ω nk s + ω 2 nk = (s + λ k1 )(s + λ k2 ) = 0 has negative real roots λ k1 and λ k2 , then (14) yields an overdamped or a critically damped behavior.
The model (10) leads to a single trajectory, while (14) always gives a continuous trajectory passing through any initial condition. Therefore, the reaching period is avoided and chattering is attenuated. Indeed, when perturbations arise, the actual robot trajectory moves away from the desired one. If a first-order sliding model is used, then the controller reacts, providing opposite phase plane trajectories towards the desired one. As some delay is inherent to the digital control, a "switching" between the curves arises, originating chattering. Using second-order curves, there is always a trajectory containing a given state, and after a perturbation, the system is not forced to follow the initial trajectory. Instead, it will follow a new one that contains the present state. As a result, the controller uses a new curve, almost parallel with the previous one, passing through the actual phase plane point. To accomplish this, the algorithm requires a second-order derivative, that is to say, it requires acceleration sensors (for more details, see [14]).
To avoid this problem, introducing an integral action was also proposed, giving rise to the reference [67]: and the Laplace expression , with two zeros and one pole at s = 0. In [68], the problem was rewritten in the scope of FC with the reference model (10) formulated as: The characteristic polynomial is now (s + λ k ) s α + λ α k , and the fractional order of −1 ≤ α ≤ 1 was explored.
Equation (16) is an "interpolation" between (14) and (15). We verified that: (i) for each robot link, Expression (16) has two zeros or two zeros and one pole, for α > 0 or α < 0, respectively; (ii) the dynamics of the manipulator, from torque to position, presents two poles; (iii) coupling phenomena can be regarded as perturbations, p k (see Figure 1). The location of the robot poles varies considerably [69], and the VSC adapts the gain to guaranty that the system's global dynamics follows the desired reference. Figure 2 sketches the root locus of the unit-feedback VSC-controlled robot links for various values of α [70,71]. The zeros and poles of the reference model must comply with the nature of the system dynamics. Expression (16) along with the manipulator dynamics establishes a compromise between zeros and poles, yielding a perceptive reference model. However, we verified that (16) requires either numerical fractional-order differentiation or fractional-order sensors.

Fractional-Order Sensor
In the classical form, an accelerometer measures the second time derivative of the displacement of a rigid body. The sensor includes a moving mass that is interconnected with a casing by one spring and one damper. Therefore, the motion of the mass with respect to the casing is able to capture the acceleration of the body.
Modern accelerometers adopt MEMSs. These sensors synthesize mechanical and electrical components in a small-scale semiconductor. MEMS technology allows integrating the three mechanical elements of the device with displacement sensors and electronics.
A uni-axial accelerometer is modeled as: where M, B, and K represent the mass, damper, and spring. The motion of M with respect to the casing is denoted by y and the absolute displacement of the casing by u.
Applying the Laplace transform, we obtain the transfer function: Equation (17) can be generalized to a fractional order as: yielding the transfer function: For γ = 2, λ = 1, and ν = 0, Equations (19) and (20) yield the dynamics of the classical accelerometer. However, (19) and (20) have no known physical meaning, and their implementation requires fractional elements that are presently not available or are unfeasible [72]. Indeed, (19) and (20) imply the fractionalization of Newton's second law, and its meaning is somewhat controversial. On the one hand, it is straightforward to generalize (17) and (18) from integer to fractional if we adopt an abstract perspective. However, on the other hand, we do not have a guarantee of the feasibility and, furthermore, that the physical properties remain of the same type. For instance, while for the standard mass, we have F(s) = Ms 2 U(s), with F(s) denoting force, meaning it has the property of being undeformable, there is no guarantee that the fractional mass, with model F(s) = Ms γ U(s), has the same behavior.
In [45], a modular N-stage cascade fractional sensor was proposed. Each k-th stage, k = 1, · · · , N, comprises the standard elements {mass, spring, damper}, represented by {M k , K k , B k }. The displacement of M k relative to the previous stage k − 1 is denoted by y k , while the external casing of the sensor corresponds to k = 0. The variables u and z k = u − ∑ k p=1 y p stand for the displacements of the device casing and M k , respectively, with respect to the external inertial frame.
The system model, having u as the input and y 1 as the output, is expressed in the Fourier domain as: where ω = 2π f is the angular frequency ( f is the frequency) and Z(jω) is given by: with Z k (jω) = jω jωB k +K k and Y k (jω) = jωM k . For the case of selecting the parameters recursively, M k+1 = η M k , B k+1 = B k and K k+1 = κK k , η, , κ ∈ R + , we obtain [43,[73][74][75]: The interval for which G N (jω) exhibits fractional dynamics is established by the relationships |jω| < B 1 M 1 and |jω| . We note that for η > 1, > 1, and η κ > 1, the influence of K k for the fractional order α is minimal [43,74,75]. Nonetheless, the recursive dependence revealed problems for small α, with arg {G N (jω)} oscillating. This problem limits the fractional behavior to the range 0 ≤ α ≤ 1, and an alternative computational procedure was proposed. Indeed, the determination of the values of the elements {M k , B k , K k }, k = 1, · · · , N, can be viewed as an optimization problem [45] with objective function: It should be noted that in this algorithm, we did not choose the fractional order, nor {ω 1 , ω 2 } in advance. The numerical scheme involves running the optimization procedure a number of times and characterizing each solution through the indices {I, σ, Ω}, where I is the fitness function value, σ stands for the standard deviation of arg {G N (jω)} in the interval ω, and Ω = log ω 1 ω 2 . The user has to choose the solution taking into account the indices. Figure 3 illustrates the frequency response of a six-stage fractional sensor. The results correspond to two instances synthesized by a PSO algorithm with a population of twenty elements and a number of iterations equal to five, yielding the fractional orders α = {0.16, 0.40}. The initialization was random in the interval {M k , B k , K k } ∈ [0, 1]. The values of the sensor parameters are listed in Table 1.

The Multidimensional Scaling Technique
MDS is a numerical procedure adopted to reduce the dimensionality and visualize high-dimensional datasets.
Let v i , i = 1, . . . , N, be objects in a space with L dimensions. First, we choose a distance d(v i , v j ), i, j = 1, . . . , N, between the pairs of objects i and j, and calculate a dissimilarity matrix D = [d(v i , v j )]. Then, we feed the MDS with D. The algorithm finds the coordinates of the objects, z i , in an embedding P-dimensional space (P ≤ L) that minimize a fitness function. The result is a matrix Z = [d(z i , z j )] that approximates D. Often, the stress cost function, S, is used as the fitness function: Nevertheless, other criteria are possible, such as the Sammon function: The MDS results were evaluated by comparing the object representations in the original and the embedding spaces. This involved the generation of the Shepard diagram, which relates d(v i , v j ) versusd(z i , z j ), and the stress map, which depicts S versus P.
A number of distances d(v i , v j ) are possible to construct D [76]. Here, we adopted {d 1 , d 2 } to quantify the dissimilarities in pairs (i, j) of objects that possess real and imaginary components. As such, the i-th object is represented by a matrix of dimension L × 2, v i = Re{v i1 }, · · · , Re{v iL } T , [Im{v i1 }, · · · , Im{v iL }] T , where Re{·} and Im{·} represent the real and imaginary parts. The distances {d 1 , d 2 } are given by: It should be noted that using (28) and (29), the objects (i, j) can be compared in the time or in the frequency domain. In the time domain, the vectors v i and v j have real components, and the distances {d 1 , d 2 } yield the standard {arccosine, Jaccard} [76]. In the frequency domain, the vectors v i and v j have complex components, and Expressions (28) and (29) are calculated directly.
In the experiments, we considered the manipulator moving from an initial state given by [7,8]: to a final one described by: To approximate the fractional derivative, we used the Al-Alaoui operator and T = 10 −4 s. In the reference model, we considered λ k = 10 (k = 1, 2), and in the control action, we adopted the saturation torque values D 1 = 200, D 2 = 100.
For each test case, we gathered the sampled state θ(t r ) = [θ 1 (t r ), θ 2 (t r )] andθ(t r ) = [θ 1 (t r ),θ 2 (t r )], t r ∈ [0, 20] s, which yields state vectors of 20,001 dimensions. The data were then organized in an N × L = (20 × 20) × (4 × 20,001)-dimensional array: We normalized the array W (t) by the average and standard deviation, µ(·) and σ(·), respectively, so that numerical saturation was avoided. Therefore, the columns of the matrix W (t) , denoted by u (t) l , with l = 1, . . . , L, are recalculated as: and a normalized array W (t) is obtained. Afterwards, the rows of W (t) , denoted by v (t) i , with i = 1, . . . , N, are used to compute the dissimilarity matrices j )], i, j = 1, . . . , 400, with h = 1, 2, that feed the MDS. Figure 5 sketches the 3-dim loci generated when adopting the arccosine distances and j ). The points, which represent the test cases, are linked by means of lines of constant δ (or α) values, while each color corresponds to one value of α (or δ). When the parameters varied, we obtained patterns that characterized the dynamics of the controlled system. For constant and low α, as δ grows, the loci describe a path that reaches an inflection zone and change direction. The round trips are identical for low values of α. As α increases, the inflection tends to vanish, and the locus becomes more insensitive to δ. For δ constant, we also verified the emergence of patterns, being similar to each other for all values of α.
The points represent the test cases. On the left, lines connect points of constant δ, and each color corresponds to points of constant α. On the right, lines connect points of constant α, and each color corresponds to points of constant δ. The fractional order α q ∈ [−0.5, 1], the width of the band δ q ∈ [10 −4 , 10 1 ], and f r ∈ [10 −2 , 10 2 ] Hz. Figure 7 portrays the 3-dim loci generated with the distances j ). We verified the existence of patterns that characterized the behavior of the system as the values of the parameters varied. For constant and low α, as δ increases, the paths spread. The round trips are similar to each other for low values of α. As α increases, the inflection tends to vanish, and the locus becomes more insensitive to δ. For constant δ, patterns also emerge, and we verified that they were of the same type, independent of the values of α. Figure 8 depicts the MDS 3-dim loci using d 1 . The points represent the test cases. On the left, lines connect points of constant δ, and each color corresponds to points of constant α. On the right, lines connect points of constant α, and each color corresponds to points of constant δ. The fractional order α q ∈ [−0.2, 0.6], the width of the band δ q ∈ [10 −4 , 10 1 ], and f r ∈ [10 −2 , 10 2 ] Hz.
To sum up, the MDS technique is able to represent large sets of objects with different classes of performance indices. The interpretation of the loci is based on the emerging patterns and clusters. The goal herein was to investigate the use of MDS to assess and easily visualize the performance of fractional-order VSC applied to revolute planar robotic manipulators. Indeed, classical indices can be adopted to assess the performance. For instance, standard time-domain parameters, such as the rise, peak and settling time, overshoot, and peak value have been widely used. However, often, single (or small set) criteria are insufficient, and different indices can lead to conflicting information. Therefore, assessing control systems' performance is a multidimensional problem that can take advantage of present-day computational information visualization techniques. The analysis adopted was based on data from the control system operation. The proposed MDS approach assumes that all relevant factors are implicitly embedded in these data. Obviously, the MDS analysis and visualization are not independent of the controlled system. However, this is not a drawback of the method. On the contrary, other systems and the influence of distinct parameters can be easily assessed by generating MDS loci that embed their corresponding information. In conclusion, the MDS constitutes a computational tool capable of generating useful representations of large sets of objects without being limited to a single class of performance indices. In this paper, the objects were closed-loop controlled systems, which differed from each other by a number of parameters. The relationships between such objects were unveiled based on the clusters and patterns that emerged on the MDS loci. More parameters can be assessed easily, generating additional objects for visualization. Therefore, further research using distinct indices and algorithms can be performed. Therefore, the results point toward further advances using other indices and fractional algorithms.

Conclusions
We addressed the use of MDS to study fractional-order VSCs. The approach compares patterns emerging in the MDS loci when varying the controller parameters. Different indices for processing the system state in the time and the frequency domains were adopted. The results illustrated the effectiveness of the MDS representation and revealed that it is simple to generate and interpret different controlled systems for a variety of test conditions. The generalization to other systems and the adoption of new assessment indices is straightforward.