A stable proportional–proportional integral tracking controller with self-organizing fuzzy-tuned gains for parallel robots

Abstract

Parallel robots are nowadays used in many high-precision tasks. The dynamics of parallel robots is naturally more complex than the dynamics of serial robots, due to their kinematic structure composed by closed chains. In addition, their current high-precision applications demand the innovation of more effective and robust motion controllers. This has motivated researchers to propose novel and more robust controllers that can perform the motion control tasks of these manipulators. In this article, a two-loop proportional–proportional integral controller for trajectory tracking control of parallel robots is proposed. In the proposed scheme, the gains of the proportional integral control loop are constant, while the gains of the proportional control loop are online tuned by a novel self-organizing fuzzy algorithm. This algorithm generates a performance index of the overall controller based on the past and the current tracking error. Such a performance index is then used to modify some parameters of fuzzy membership functions, which are part of a fuzzy inference engine. This fuzzy engine receives, in turn, the tracking error as input and produces an increment (positive or negative) to the current gain. The stability analysis of the closed-loop system of the proposed controller applied to the model of a parallel manipulator is carried on, which results in the uniform ultimate boundedness of the solutions of the closed-loop system. Moreover, the stability analysis developed for proportional–proportional integral variable gains schemes is valid not only when using a self-organizing fuzzy algorithm for gain-tuning but also with other gain-tuning algorithms, only providing that the produced gains meet the criterion for boundedness of the solutions. Furthermore, the superior performance of the proposed controller is validated by numerical simulations of its application to the model of a planar three-degree-of-freedom parallel robot. The results of numerical simulations of a proportional integral derivative controller and a fuzzy-tuned proportional derivative controller applied to the model of the robot are also obtained for comparison purposes.

Introduction

Parallel robots are nowadays used in many high-precision industrial tasks, such as machining, ^1,2 welding, ^3,4 packaging, ^5,6 as well as flight simulators ⁷ and telescopes. ^8,9 These high-precision applications demand the implementation of position regulators, for point-to-point motion tasks, and tracking controllers, when the end-effector motion has to follow a prescribed trajectory of motion. ¹⁰ Moreover, parallel robots consist of closed kinematic chains that can be modeled, roughly speaking, as a collection of open kinematic chains in which their extremes (at the end-effectors) are joined together. ¹¹ This parallel structure can be mathematically represented as a dynamic model of a serial robot plus a set of holonomic constraints. Therefore, as it will be explained later, the dynamic model of parallel robots has inherently a more complex structure than the dynamic model of serial robots. This has led researchers in robot control to conceive more robust controllers to cope with inaccuracies and parameter uncertainties that may exist in the dynamic models of parallel robots. Hence, we can see that the development of trajectory tracking controllers for parallel robots is nowadays an important and interesting problem.

Recent approaches for control of parallel manipulators found in the literature include model-based controllers ^12,13 and not model-based controllers. ^12,14,15 In the study by Ren et al., ¹² a comparison of several control approaches for robot tracking of a three-degree-of-freedom (3-DOF) parallel robot is presented. They compare the performance of an adaptive controller, a proportional integral (PI)-type synchronized controller (model-based), a conventional proportional integral derivative (PID) controller, and an adaptive synchronized controller (not model-based controller). In Diaz-Rodriguez et al.’s ¹⁶ work, a reduced model-based control of a 3-DOF parallel robot is proposed. The reduced model is obtained by considering a simplified model with a set of relevant parameters. In the study by Borbonnais et al., ¹³ a computer torque controller and a conventional PID controller are implemented for a novel five-bar parallel robot.

An interesting two-loop control scheme is proposed by Kelly and Moreno ¹⁷ for motion control of manipulators in operational space. The proposed scheme is composed of a joint velocity inner loop and an operational space kinematic control outer loop. The motivation of the cited work is the fact that in many industrial robots the motion control of each joint is performed using inner joint velocity control loops and outer position control loops. ^{17 –19} In addition, it is shown that the overall closed-loop dynamics of the two-loop control scheme can be seen as a cascade system, and it presents sufficient conditions for the Global Exponential Stability of the closed-loop systems. Inspired by Kelly and Moreno’s work, ¹⁷ Soto and Campa ²⁰ proposed a two-loop control scheme for redundant manipulators. In such a scheme, an outer kinematic controller generates the joint velocity commands for an inner velocity controller for each joint of the manipulator. A recent application of a two-loop motion control scheme for a parallel robot is reported by Campa et al., ²¹ where a kinematic model-based operational space tracking controller with a two-loop structure for a hexapod-type parallel robot is proposed. The two-loop structure contains a Resolved Motion Rate Controller (RMRC) in the outer loop and a joint velocity PI in the inner loop.

On the other hand, not model-based controllers are not sensitive to model uncertainties or inaccuracies. This can be an important advantage when dealing with the dynamic model of parallel robots, which have inherently a more complex structure than the dynamic model of serial robots, as it was pointed above. Regarding classic not model-based controllers, such as PD or PID, these can improve their robustness when used in variable gains schemes or when a nonlinear function of the error is added. In these schemes, either the controller gains or some nonlinear control action are set-up online according to a measure of the performance of the controller. In the study by Salinas et al., ¹⁵ a family of nonlinear PID-like controllers is proposed in which an integral action of a nonlinear function of the position error is added to the control signal. Recent approaches of variable gains control schemes for parallel robots include ¹⁴ in which a discrete-time motion controller is proposed. In such work, the output of an Adaptive Neuro-Fuzzy Inference System is used in the control law in addition to a conventional PID controller. In Villareal-Cervantes and Alvarez-Gallegos’ work, ²² a PID control off-line tuning scheme using Evolutionary Algorithms is proposed. The scheme is tested on a planar parallel robot with a five-bar mechanism. In Ouyang et al.’s work, ²³ a nonlinear PD control for a parallel manipulator system is proposed, in which the P and D gains are certain class of nonlinear bounded functions of the position error. In the study by Zhu and Zhang, ²⁴ a fuzzy-tuned PID control strategy is proposed for regulation control of a parallel robot. In Stan et al.’s work, ²⁵ a fuzzy controller performs the fine tuning of the gains of a PID controller. A coarse tuning of these gains is previously done. The fuzzy controller is applied to a Cartesian parallel robot of 3-DOF known as ISOGLIDE. In another related recent work, ²⁶ a fuzzy-PID controller is designed for control of a 6-DOF hydraulically driven parallel robot. The PID gains are tuned by a conventional fuzzy algorithm which uses the joint position errors and its temporal derivatives as input variables. However, stability analysis of these variable gains schemes for not model-based controllers for parallel robots have not been reported in the literature. The stability of the resultant closed-loop system is an important concern in the design process of controllers for manipulators.

Given the larger complexity of the dynamic models of parallel robots, and the need of improved robustness of the controllers before parameter variations or uncertainties, in this work, we propose a novel design of a two-loop P-PI controller for tracking control of parallel manipulators. The outer loop is a joint position P control loop, and the inner loop is a joint velocity PI control loop. The gains of the inner PI control loop are constant, while the gains of the outer P control loop are tuned by a self-organizing fuzzy (SOF) algorithm, which is a novel algorithm that generates the gain depending on the current tracking errors. This algorithm employs a set of fuzzy rules that can be updated online according to a performance index of the controller response. This gain tuning scheme provides the needed increased robustness in presence of perturbations and parameter uncertainty of the overall controller, as well as improved results when it is tracking fast changing trajectories. The usefulness and performance of the proposed controller is tested by numerical simulations, in comparison with two controllers, a PID controller, and a fuzzy-tuned PD (FPD) controller. Moreover, the stability analysis of the closed-loop system is carried out, giving as a result the proof of the Uniform Ultimate Boundedness of the solutions of the closed-loop system.

The main contributions of this work are the design and the stability analysis of a tracking control structure based on a two-loop scheme composed of two not model-based controllers, that is, P and PI. In addition, the gains of the P control loop can be tuned online. Moreover, the P-gain tuning is performed by a novel SOF algorithm. Furthermore, the stability analysis of the resulting closed-loop system is carried out. It is worth to mention that the result of the stability analysis holds for any tuning algorithm for tuning the P gains, by only ensuring that some conditions on the boundedness of the produced gains are satisfied. Therefore, this is a general result for a family of two-loop P-PI tracking controllers in which the P gain of the outer control loop is variable. According to the best knowledge of the authors, this is an original result not published before elsewhere.

The remaining of this article is organized as follows. In the “Preliminaries” section, the dynamic model of the parallel robot as well as some of their properties useful for the stability analysis are described. In the section “Proposed controller,” the full description of the proposed controller is given. In section “Stability analysis,” the development of the stability analysis in the sense of Lyapunov of the closed-loop system is carried out, and the result of this analysis is presented. In section “Simulations,” the elements of the dynamic and kinematic model of a 3-DOF parallel robot are described as well as the desired trajectory for simulations. Also, the numerical simulations conducted are presented and discussed. Finally, some concluding remarks are given.

Preliminaries

Notation

In this work, vectors are denoted with italic-bold lowercase letters, for example, $\mathit{x}$ or $\mathit{\omega }$ . Matrices are denoted with capital letters, for example, A. $\parallel \mathit{x}\parallel =\sqrt{{\mathit{x}}^{\text{T}}\mathit{x}}$ represents the Euclidean norm of vector $\mathit{x}$ . ${\lambda }_{\text{max}}\left\{A\right\}$ and ${\lambda }_{\text{min}}\left\{A\right\}$ represent the largest and smallest eigenvalues of matrix A, respectively.

Dynamic model of the robot

The joints of a parallel manipulator can be classified as n actuated joints, represented by vector $\mathit{q}\in {ℝ}^n$ , and m non-actuated joints, represented by vector $\beta \in {ℝ}^m$ . Moreover, the l operational or Cartesian variables of the end-effector can be represented with vector $\mathit{x}\in {ℝ}^l$ . Let us define the vector of constrained variables $\rho$ as

$\mathit{\rho }={\left[{\mathit{q}}^T{\beta }^T{\mathit{x}}^T\right]}^{\text{T}}\in {ℝ}^s$

with $s=n+m+l$ . By applying the Euler–Lagrange formulation, the dynamic model of a parallel robot ²⁷ in absence of friction is in general formulated as

$M^{\prime }(\mathit{\rho })\ddot{\mathit{\rho }}+C^{\prime }(\mathit{\rho },\dot{\mathit{\rho }})\dot{\mathit{\rho }}+{\mathit{g}}^{\prime }(\mathit{\rho })={\mathit{\tau }}_{\mathit{\rho }}+D^{\text{T}}(\mathit{\rho })\mathit{\lambda }$

1

where $M^{\prime }(\mathit{\rho })\in {ℝ}^{s\times s}$ is the symmetric positive definite inertia matrix, $C^{\prime }(\mathit{\rho },\dot{\mathit{\rho }})\in {ℝ}^{s\times s}$ is the Coriolis and centripetal forces matrix, ${\mathit{g}}^{\prime }(\mathit{\rho })\in {ℝ}^s$ is the vector of gravitational forces, ${\mathit{\tau }}_{\mathit{\rho }}\in {ℝ}^s$ is the vector of generalized forces associated to scalar variables of $\mathit{\rho }$ , $D(\mathit{\rho })=\partial \mathit{\gamma }(\mathit{\rho })/\partial \mathit{\rho }\in {ℝ}^{(k)(n)\times (s)}$ is the Jacobian matrix of the system holonomic constraints $\mathit{\gamma }(\mathit{\rho })\in {ℝ}^{(k)(n)}$ , with k as the minimum number of coordinates in the space that describes the position of the end-effector, and $\mathit{\lambda }\in {ℝ}^{(k)(n)}$ is the vector of Lagrange multipliers.

Kinematic transformation

Let us define a kinematic transformation between vectors $\dot{\mathit{\rho }}$ and $\dot{\mathit{q}}$ as

$\dot{\mathit{\rho }}=R(\mathit{\rho })\dot{\mathit{q}}$

2

to obtain a so-called reduced dynamic model. Notice that the matrix $R(\mathit{\rho })\in {ℝ}^{(s)\times (n)}$ can be constructed with the Jacobian matrices $J_{\beta }(\mathit{\rho })$ and $J_{\mathit{x}}(\mathit{\rho })$ from the differential kinematic models $\dot{\beta }=J_{\beta }(\mathit{\rho })\dot{\mathit{q}}$ and $\dot{\mathit{x}}=J_{\mathit{x}}(\mathit{\rho })\dot{\mathit{q}}$ , as

$R(\mathit{\rho })=\left[\begin{array}{c}I\\ J_{\beta }(\rho )\\ J_{\mathit{x}}(\rho )\end{array}\right],I\in {ℝ}^{n\times n}$

The derivative with respect to time of equation (2) is

$\ddot{\mathit{\rho }}=\dot{R}(\mathit{\rho })\dot{\mathit{q}}+R(\mathit{\rho })\ddot{\mathit{q}}$

3

By substituting $\dot{\mathit{\rho }}$ from equation (2) and $\ddot{\mathit{\rho }}$ from equation (3) in equation (1), we get

$M^{\prime }(\mathit{\rho })\dot{R}(\mathit{\rho })\dot{\mathit{q}}+M^{\prime }(\mathit{\rho })R(\mathit{\rho })\ddot{\mathit{q}}+C^{\prime }(\mathit{\rho },\dot{\mathit{\rho }})R(\mathit{\rho })\dot{\mathit{q}}+{\mathit{g}}^{\prime }(\mathit{\rho })={\mathit{\tau }}_{\mathit{\rho }}+D^{\text{T}}(\mathit{\rho })\mathit{\lambda }$

4

After multiplying equation (4) by $R^T(\dot{\mathit{\rho }})$ and rearranging terms we obtain

$\begin{array}{l}{R}^{\text{T}}(\mathit{\rho })M^{\prime }(\mathit{\rho })R(\mathit{\rho })\ddot{\mathit{q}}\\ +\left[R^{\text{T}}(\mathit{\rho })M^{\prime }(\mathit{\rho })\dot{R}(\mathit{\rho })+R^{\text{T}}(\mathit{\rho })C^{\prime }(\mathit{\rho },\dot{\mathit{\rho }})R(\mathit{\rho })\right]\dot{\mathit{q}}\\ +R^{\text{T}}(\mathit{\rho }){\mathit{g}}^{\prime }(\mathit{\rho })=R^{\text{T}}(\mathit{\rho }){\mathit{\tau }}_{\mathit{\rho }}\end{array}$

5

where the term with the product of the Jacobian matrix of the constraints by the Lagrange multipliers vanishes because it belongs to the null space of $D^{\text{T}}(\mathit{\rho })$ . By defining

$\begin{array}{c}M(\mathit{\rho })=R^{\text{T}}(\mathit{\rho })M^{\prime }(\mathit{\rho })R(\mathit{\rho })\\ C(\mathit{\rho },\dot{\mathit{\rho }})=R^{\text{T}}(\mathit{\rho })M^{\prime }(\mathit{\rho })\dot{R}(\mathit{\rho })+R^{\text{T}}(\mathit{\rho })C^{\prime }(\mathit{\rho },\dot{\mathit{\rho }})R(\mathit{\rho })\\ \mathit{g}(\mathit{\rho })=R^{\text{T}}(\mathit{\rho }){\mathit{g}}^{\prime }(\mathit{\rho })\\ \mathit{\tau }=R^{\text{T}}(\mathit{\rho }){\mathit{\tau }}_{\mathit{\rho }}\end{array}$

equation (5) can be rewritten as

$M(\mathit{\rho })\ddot{\mathit{q}}+C(\mathit{\rho },\dot{\mathit{\rho }})\dot{\mathit{q}}+\mathit{g}(\mathit{\rho })=\mathit{\tau }$

6

In Ghorbel et al.’s work, ¹¹ it is proven, by means of the implicit function theorem, that there exists a unique parametrization of $\mathit{\rho }\in {\mathcal{N}}_{\mathit{\rho }}$ , inside a neighborhood ${\mathcal{N}}_{\mathit{\rho }}$

$\mathit{\rho }=\eta (\mathit{q})$

for any $\mathit{q}\in {\mathcal{N}}_{\mathit{q}}$ inside a neighborhood ${\mathcal{N}}_{\mathit{q}}$ , whenever the system is not in a singular configuration.

In addition, Muller ²⁸ established that, for a parallel machine, a subset $\mathit{q}$ of δ joint variables determines its configuration, in virtue of that there exists a smooth mapping ψ that assigns to each $\mathit{\rho }$ , the parallel machine configuration as $\mathit{\rho }=\psi (\mathit{q})$ , where the map ${\psi }^{-1}$ is a local parametrization of the δ dimensional manifold $\mathcal{V}$

$\mathcal{V}=\left\{\mathit{\rho }\in {\mathcal{V}}^n;\mathit{\gamma }(\mathit{\rho })=0\right\}$

where $\mathcal{V}$ represents the set of all admissible configurations of the parallel machine, and $\mathit{\gamma }(\mathit{\rho })$ represents the holonomic constraints.

In consequence, without loss of generality, we can write down the matrices of the dynamic model $M(\mathit{\rho }),C(\mathit{\rho },\dot{\mathit{\rho }})$ , and $\mathit{g}(\mathit{\rho })$ as $M(\mathit{q})$ , $C(\mathit{q},\dot{\mathit{q}})$ , and $\mathit{g}(\mathit{q})$ , respectively. Thus, the dynamic model in equation (6) takes the form

$M(\mathit{q})\ddot{\mathit{q}}+C(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{g}(\mathit{q})=\mathit{\tau }$

7

where $M(\mathit{q}),C(\mathit{q},\dot{\mathit{q}})$ , and $\mathit{g}(\mathit{q})$ are the inertia matrix, the Coriolis and centripetal forces matrix, and the gravitational forces vector of the reduced model, respectively, and $\mathit{\tau }$ is the vector of forces associated with the actuated joints $\mathit{q}$ . The model in equation (7) exhibits the following properties.

Property 1

M(q) is symmetric and positive definite. ²⁹

Property 2

$\dot{M}(\mathit{q})-2C(\mathit{q},\dot{\mathit{q}})$ is skew-symmetric. ^11,29

Property 3

There exists a constant $k_C>0$ such that $\parallel C^{\prime }(\mathit{q},\dot{\mathit{q}})\parallel \le k_C\parallel \dot{\mathit{q}}\parallel$ , for all $\mathit{q}\in {ℝ}^n$ . ³⁰

Property 4

For robots with only revolute joints, there exists a constant $k^{\prime }$ such that $\parallel \mathit{g}(\mathit{q})\parallel \le k^{\prime }$ , for all $\mathit{q}\in \Omega \subset {ℝ}^n$ .

Proof

For robots with only revolute joints, the entries of the vector of gravitational forces ${\mathit{g}}^{\prime }(\mathit{\rho })$ (see equation (1)) are sinusoidal or trigonometric functions with finite coefficients. Thus the entries are upper bounded. Now, lets consider the vector $\mathit{g}(\mathit{\rho })=R^{\text{T}}(\mathit{\rho }){\mathit{g}}^{\prime }(\mathit{\rho })$ . Note that the entries of this vector are upper bounded for all $\mathit{q}\in \Omega \subset {ℝ}^n$ , where Ω denotes a space of configurations of the robot in which the matrix $R(\mathit{\rho })$ is well posed.

Proposed controller

In this work, we propose a controller for trajectory tracking of parallel manipulators. The structure of the proposed controller is a two-loop cascade structure, with a P position controller in the outer loop and a PI controller in the inner loop. The overall structure of the controller is shown in Figure 1. Note that the desired joint position and velocity vectors ${\mathit{q}}_d$ and ${\dot{\mathit{q}}}_d$ , respectively, are represented in the Figure 3, as well as the position error and velocity error vectors $\tilde{\mathit{q}}={\mathit{q}}_d-\mathit{q}$ and $\dot{\tilde{\mathit{q}}}={\dot{\mathit{q}}}_d-\dot{\mathit{q}}$ , respectively. We define the velocity control law of the outer loop as

${\mathit{\omega }}_d=\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}+{\dot{\mathit{q}}}_d$

8

Figure 1.

Proposed controller. SOF: self-organizing fuzzy.

Notice that $\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})$ is a diagonal matrix of variable gains of $n\times n$ dimension. The values of each diagonal element of the matrix are tuned by an algorithm, represented in Figure 1 as the block labeled SOF controller, that evaluates the performance of the current gain values by measuring the position and velocity errors. Such an algorithm will be described in a later section in this article. Regarding the Figure 1, the velocity error $\tilde{\mathit{\omega }}$ in the inner loop can be written as

$\begin{array}{c}\tilde{\mathit{\omega }}=\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}+{\dot{\mathit{q}}}_d-\dot{\mathit{q}}\\ =\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}+\dot{\tilde{\mathit{q}}}\end{array}$

9

where ${\dot{\mathit{q}}}_d-\dot{\mathit{q}}$ has been substituted by $\dot{\tilde{\mathit{q}}}$ . Notice that, in Figure 1, the diagonal arrow shown over the block labeled as $\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})$ represents the tuning of the gains by the SOF algorithm.

In the inner loop, the control law can be written as

$\mathit{\tau }=K_{vp}\tilde{\mathit{\omega }}+K_{vi}\xi$

10

where

$\dot{\xi }=\tilde{\mathit{\omega }},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\xi ={\displaystyle {\int }_0^t}\hspace{0.17em}\tilde{\mathit{\omega }}(\mathit{\sigma })\hspace{0.17em}\text{d}\mathit{\sigma }$

11

$K_{vp}$ and $K_{vi}$ are $n\times n$ diagonal matrices of constant gains. These matrices can be defined for a later analysis as

$K_{vp}:=k_{p}I$

12

$K_{vi}:=\alpha K_{vp}$

13

where $k_p$ and α are positive constants and $I\in {ℝ}^n$ is the identity matrix.

SOF controller

The proposed algorithm for tuning the gains of the outer P loop is a novel SOF controller. This controller generates the gains based on the current position and velocity errors. The performance of the overall control strategy is evaluated to produce a performance index. This parameter is then used to adjust the parameters of the membership functions of the system to obtain an improved control signal. In a fuzzy controller, the membership functions are functions that map a crisp input variable to a membership grade that can take continuous values in the interval $\left[0,\dots \mathrm{,1}\right]$ . ³¹ This process is known as fuzzification. In summary, by continuously evaluating the performance of the controller and generating the respective changes in the fuzzy membership parameters, the performance of the tracking is improved. A simplified block diagram of the structure of the SOF controller is shown in Figure 2. The input signals are the ith position and velocity errors ${\tilde{\mathit{q}}}_i$ and ${\dot{\tilde{\mathit{q}}}}_i$ . The output is ${\varphi }_i$ , which is the $i\text{th}$ diagonal element of variable gains matrix Φ. In the following subsections, the stages of the SOF controller are described.

Figure 2.

Simplified structure of the SOF controller. SOF: self-organizing fuzzy.

Input stage

A block diagram of the structure of the input stage is shown in Figure 2. In this stage, two processes are performed: scaling and quantization, to normalize the levels of input signals before passing them to the following stage. The scaling is done by multiplying by a scaling factor. Next, the scaled signal is subject to a quantization process by mapping it into two sets of discrete values $S_R=\left\{1,2,3,4\right\}$ and $S_F=\left\{-6,-5.5,-5,\dots \mathrm{,5,5.5,6}\right\}$ . These two discrete sets correspond to two types of discretization which are conveniently labeled regular and fine, respectively. Moreover, the values of the discrete set $S_R$ corresponds to the linguistic variables zero, small, medium, and large, labeled with letters Z, S, M, and L. Position and velocity errors ${\tilde{\mathit{q}}}_i$ and ${\dot{\tilde{\mathit{q}}}}_i$ are scaled and quantized. In addition, the output ${\varphi }_i$ is also fed back into the scaling and quantization processes. The output signals from regular quantization of ${\tilde{\mathit{q}}}_i$ and ${\dot{\tilde{\mathit{q}}}}_i$ , labeled as $e_{p_I}$ and $e_{v_I}$ , respectively, are used to address a performance index table. This table, which is described in a subsequent section, produces an appropriate value of the performance index.

On the other hand, the output signals from fine quantization of ${\tilde{\mathit{q}}}_i$ and ${\dot{\tilde{\mathit{q}}}}_i$ , labeled as $e_{p_F}$ and $e_{v_F}$ , respectively, are submitted to a fuzzification process. The output signal from fine quantization of ${\varphi }_i$ is stored in a First Input–First Output (FIFO) serial-shifting memory. In this memory, each introduced value is stored during four time periods (see Figure 3). Then, the sum of the stored value and the current performance index is stored into a four-element serial-shifting memory. The output of these four memory elements are used as the centers of four fuzzy membership functions.

Figure 3.

Structure of the input stage.

Fuzzification

The fuzzy membership functions selected are triangular shaped functions, with constant width and adjustable central values. These functions are shown in Figure 4. These central values are considered as the inference rules of the fuzzy mechanism. When crisp values $e_{p_f}$ and $e_{v_f}$ are fuzzified by the four membership functions, the resulting memberships grades are the sets $\left\{{\mu }_{1p},{\mu }_{2p},{\mu }_{3p},{\mu }_{4p}\right\}$ and $\left\{{\mu }_{1v},{\mu }_{2v},{\mu }_{3v},{\mu }_{4v}\right\}$ , respectively. Next, these membership grades are fed into a fuzzy inference process that uses the min function for each pair of fuzzy membership grades of position and velocity errors. The output $I_K$ from this inference algorithm is the vector

Figure 4.

Fuzzification of discretized errors.

$I_K=[\text{min}({\mu }_{1p},{\mu }_{1v}),\text{min}({\mu }_{2p},{\mu }_{2v}),\text{min}({\mu }_{3p},{\mu }_{3v}),\text{min}({\mu }_{4p},{\mu }_{4v})]$

Defuzzification

Defuzzification is carried out by the Mean of Maxima algorithm. In this algorithm, the two greatest values of vector $I_K$ are multiplied by the central value of the fuzzy membership function from which the membership grade is obtained. The output $U_i$ is the arithmetic average of both products. Finally, this output is descaled by multiplying by a reverse scaling gain $k_r$ , and the product is added to the past period gain

${\varphi }_i(t)={\varphi }_i(t-T)+k_r{U}_i$

where ${\varphi }_i(t-T)$ is the new value of the gain and ${\varphi }_i(t-T)$ is the past value.

Performance index table

Performance index values are generated according to the current position and velocity error values. In Table 1, the selected values for the performance index are shown. These values were selected according to the following criterion: if $e_{p_I}$ or $e_{v_I}$ are Z (Zero), we can assume that the errors are decreasing. Then, the performance is considered good and the performance index is set to zero. If $e_{p_I}$ is L (Large) and $e_{v_I}$ is S (Small), we can assume, similarly, that the errors are decreasing and the performance index is set to zero. If $e_{p_I}$ and $e_{v_I}$ are S, M (Medium) or L, it is considered that the performance is poor, then the performance index is set to an appropriate positive or negative value to modify the rule, to adjust the produced gain value accordingly.

Table 1.

Performance index table.

$e_{p_I}/e_{v_I}$	Z	S	M	L
Z	0	0	0	0
S	0	−2	−2	−2
M	0	4	4	4
L	0	0	2	2

Z: zero; S: small; M: medium; L: large.

Stability analysis

In this work, the approach that we address to carry out the stability analysis of the closed-loop system is considering the SOF-P-PI controller simply as a variable gains P-PI controller, without taking into account the dynamics of the tuning algorithm used for tuning the variable gains. We only assume that the values of the tuned gains of the proportional outer loop can be upper and lower bounded, as well as the change ratio of these gains.

Closed-loop system

By substituting equation (10) in equation (7) we obtain

$M(\mathit{q})\hspace{0.17em}\ddot{\mathit{q}}+C(\mathit{q},\dot{\mathit{q}})\hspace{0.17em}\dot{\mathit{q}}+\mathit{g}(\mathit{q})=K_{vp}\tilde{\mathit{\omega }}+K_{vi}\mathit{\xi }$

Since the velocity error in the inner loop and its temporal derivative are

$\begin{array}{l}\tilde{\mathit{\omega }}={\mathit{\omega }}_d-\dot{\mathit{q}}\dot{\tilde{\mathit{\omega }}}={\dot{\mathit{\omega }}}_d-\ddot{\mathit{q}}\end{array}$

we can write

$\begin{array}{l}M(\mathit{q})\dot{\tilde{\mathit{\omega }}}+C(\mathit{q},\dot{\mathit{q}})\tilde{\mathit{\omega }}+K_{vp}\tilde{\mathit{\omega }}+K_{vi}\mathit{\xi }\\ =M(\mathit{q}){\dot{\mathit{\omega }}}_d+C(\mathit{q},\dot{\mathit{q}}){\mathit{\omega }}_d+\mathit{g}(\mathit{q})\end{array}$

14

After adding $\alpha M(\mathit{q})\tilde{\mathit{\omega }}+\alpha C(\mathit{q},\dot{\mathit{q}})\mathit{\xi }$ to both sides of equation (14), we obtain

$\begin{array}{l}M(\mathit{q})\left[\dot{\tilde{\mathit{\omega }}}+\alpha \tilde{\mathit{\omega }}\right]+C(\mathit{q},\dot{\mathit{q}})\left[\tilde{\mathit{\omega }}+\alpha \xi \right]+K_{vp}\tilde{\mathit{\omega }}+K_{vi}\mathit{\xi }\\ =M(\mathit{q})\left[{\dot{\mathit{\omega }}}_d+\alpha \tilde{\mathit{\omega }}\right]+C(\mathit{q},\dot{\mathit{q}})\left[{\mathit{\omega }}_d+\alpha \xi \right]+\mathit{g}(\mathit{q})\end{array}$

15

Let define the state vector $\mathit{z}={\left[{\tilde{\mathit{q}}}^T\hspace{0.17em}\hspace{0.17em}{\mathit{\xi }}^T\hspace{0.17em}\hspace{0.17em}{\tilde{\mathit{\omega }}}^T\right]}^T\in {ℝ}^{3n}$ . Let the term $d(t,\mathit{z})$

$\mathit{d}(t,\mathit{z})=M(\mathit{q})\left[{\dot{\mathit{\omega }}}_d+\alpha \tilde{\mathit{\omega }}\right]+C(\mathit{q},\dot{\mathit{q}})\left[{\mathit{\omega }}_d+\alpha \mathit{\xi }\right]+\mathit{g}(\mathit{q})$

16

be considered as a disturbance. From equations (9), (11), (15), and (16) the closed-loop system is

$\frac{\text{d}}{\text{d}t}\left[\begin{array}{c}\tilde{\mathit{q}}\\ \xi \\ \tilde{\mathit{\omega }}\end{array}\right]=\left[\begin{array}{c}\tilde{\mathit{\omega }}-\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}\\ \tilde{\mathit{\omega }}\\ -\alpha \tilde{\mathit{\omega }}-M^{-1}(\mathit{q})D(\mathit{q},\dot{\mathit{q}},\mathit{z})\end{array}\right]$

17

where

$D(\mathit{q},\dot{\mathit{q}},\mathit{z})=C(\mathit{q},\dot{\mathit{q}})\left[\tilde{\mathit{\omega }}+\alpha \mathit{\xi }\right]+K_{vp}\tilde{\mathit{\omega }}+K_{vi}\mathit{\xi }-\mathit{d}(t,\mathit{z})$

Note that the origin of the closed-loop system in equation (17) is not an equilibrium. In the current case, the stability of a nominal system will be investigated first, to investigate the uniform and uniform ultimate boundedness of the solutions of the closed-loop system in equation (17).

Stability of the nominal system

When the disturbance $\mathit{d}(t,\mathit{z})$ is equal to zero, the closed-loop system is a nominal system

$\frac{\text{d}}{\text{d}t}\left[\begin{array}{c}\tilde{\mathit{q}}\\ \xi \\ \tilde{\mathit{\omega }}\end{array}\right]=\left[\begin{array}{c}\tilde{\mathit{\omega }}-\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}\\ \tilde{\mathit{\omega }}\\ -\alpha \tilde{\mathit{\omega }}-M^{-1}(\mathit{q})E(\mathit{q},\dot{\mathit{q}},\mathit{z})\end{array}\right]$

18

where

$E(\mathit{q},\dot{\mathit{q}},\mathit{z})=C(\mathit{q},\dot{\mathit{q}})\left[\tilde{\mathit{\omega }}+\alpha \mathit{\xi }\right]+K_{vp}\tilde{\mathit{\omega }}+K_{vi}\mathit{\xi }$

Notice that the origin $\mathit{z}=0\in {ℝ}^{3n}$ is the only equilibrium point of equation (18). We propose the positive definite Lyapunov candidate function

$\begin{array}{l}V(\tilde{\mathit{q}},\xi ,\tilde{\mathit{\omega }})=\frac{1}{2}{\tilde{\mathit{q}}}^{\text{T}}\tilde{\mathit{q}}+\alpha {\xi }^{\text{T}}\left[k_{p}I+\frac{1}{2}\alpha M(\mathit{q})\right]\mathit{\xi }\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\alpha {\xi }^{\text{T}}M(\mathit{q})\tilde{\mathit{\omega }}+\frac{1}{2}{\tilde{\mathit{\omega }}}^{\text{T}}M(\mathit{q})\tilde{\mathit{\omega }}\end{array}$

19

The first term of equation (19) can be written as $\frac{1}{2}{\tilde{\mathit{q}}}^T\tilde{\mathit{q}}=\frac{1}{2}\parallel \tilde{\mathit{q}}{\parallel }^2$ . The second term of $V(\tilde{\mathit{q}},\mathit{\xi },\tilde{\mathit{\omega }})$ can be lower bounded as follows

$\alpha {\mathit{\xi }}^{\text{T}}\left[k_{p}I+\frac{1}{2}\alpha M(\mathit{q})\right]\mathit{\xi }\ge \left[\frac{1}{2}{\alpha }^2{\mathit{\lambda }}_{\text{min}}\left\{M\right\}+\alpha k_p\right]\parallel \mathit{\xi }{\parallel }^2$

and it can be upper bounded as

$\alpha {\mathit{\xi }}^{\text{T}}\left[k_{p}I+\frac{1}{2}\alpha M(\mathit{q})\right]\mathit{\xi }\le \left[\frac{1}{2}{\alpha }^2{\mathit{\lambda }}_{\text{max}}\left\{M\right\}+\alpha k_p\right]\parallel \mathit{\xi }{\parallel }^2$

The third and fourth terms can be lower and upper bounded as

$\begin{array}{l}-\alpha {\mathit{\lambda }}_{\text{max}}\parallel \xi \parallel \parallel \tilde{\mathit{\omega }}\parallel \le \hspace{1em}\alpha {\xi }^{T}M(\mathit{q})\tilde{\mathit{\omega }}\le \alpha {\mathit{\lambda }}_{\text{max}}\parallel \xi \parallel \parallel \tilde{\mathit{\omega }}\parallel ,\\ \frac{1}{2}{\mathit{\lambda }}_{\text{min}}\left\{M\right\}\parallel \tilde{\mathit{\omega }}{\parallel }^2\le \hspace{1em}\frac{1}{2}{\tilde{\mathit{\omega }}}^{T}M(\mathit{q})\tilde{\mathit{\omega }}\le \frac{1}{2}{\mathit{\lambda }}_{\text{max}}\left\{M\right\}\parallel \tilde{\mathit{\omega }}{\parallel }^2\end{array}$

Thus $V(\tilde{\mathit{q}},\mathit{\xi },\tilde{\mathit{\omega }})$ in equation (19) can be bounded as

${\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]}^{\text{T}}{P}_1\hspace{0.17em}\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]\le V(\tilde{\mathit{q}},\xi ,\tilde{\mathit{\omega }})\le {\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]}^{\text{T}}{P}_2\hspace{0.17em}\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]$

where

$\begin{array}{c}{P}_1=\left[\begin{array}{ccc}\frac{1}{2}& 0& 0\\ 0& \alpha k_p+\frac{1}{2}{\alpha }^2{\mathit{\lambda }}_{\text{min}}\left\{M\right\}& -\frac{1}{2}\alpha {\mathit{\lambda }}_{\text{max}}\left\{M\right\}\\ 0& -\frac{1}{2}\alpha {\mathit{\lambda }}_{\text{max}}\left\{M\right\}& \frac{1}{2}{\mathit{\lambda }}_{\text{min}}\left\{M\right\}\end{array}\right],\\ P_2=\left[\begin{array}{ccc}\frac{1}{2}& 0& 0\\ 0& \alpha k_p+\frac{1}{2}{\alpha }^2{\mathit{\lambda }}_{\text{max}}\left\{M\right\}& \frac{1}{2}\alpha {\mathit{\lambda }}_{\text{max}}\left\{M\right\}\\ 0& \frac{1}{2}\alpha {\mathit{\lambda }}_{\text{max}}\left\{M\right\}& \frac{1}{2}{\mathit{\lambda }}_{\text{max}}\left\{M\right\}\end{array}\right]\end{array}$

Then, the Lyapunov candidate function in equation (19) satisfies

${\mathit{\lambda }}_1{\Vert \begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\Vert }^2\le V(\tilde{\mathit{q}},\xi ,\tilde{\mathit{\omega }})\le {\mathit{\lambda }}_2{\Vert \begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\Vert }^2$

20

${\mathit{\lambda }}_1={\mathit{\lambda }}_{\text{min}}\left\{P_1\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\lambda }}_2={\mathit{\lambda }}_{\text{max}}\left\{P_2\right\}$

21

To ensure that the Lyapunov candidate function is positive definite, it is enough to prove that matrix $P_1$ is positive definite, which is satisfied if

$k_p<\frac{\alpha \left({\mathit{\lambda }}_{\text{max}}^2\left\{M\right\}-{\mathit{\lambda }}_{\text{min}}^2\left\{M\right\}\right)}{2{\mathit{\lambda }}_{\text{min}}\left\{M\right\}}$

22

The temporal derivative of the Lyapunov candidate function in equation (19) is

$\begin{array}{c}\dot{V}(\tilde{\mathit{q}},\xi ,\tilde{\mathit{\omega }})={\tilde{\mathit{q}}}^{\text{T}}\dot{\tilde{\mathit{q}}}+2\alpha k_p{\xi }^{\text{T}}\dot{\xi }+{\alpha }^2{\xi }^{\text{T}}M(\mathit{q})\dot{\xi }\\ +\frac{1}{2}{\alpha }^2{\xi }^{\text{T}}\dot{M}(\mathit{q})\xi +\alpha {\dot{\xi }}^{\text{T}}M(\mathit{q})\tilde{\mathit{\omega }}\\ +\alpha {\xi }^{\text{T}}\dot{M}(\mathit{q})\tilde{\mathit{\omega }}+\alpha {\xi }^{\text{T}}M(\mathit{q})\dot{\tilde{\mathit{\omega }}}\\ +{\tilde{\mathit{\omega }}}^{\text{T}}M(\mathit{q})\dot{\tilde{\mathit{\omega }}}+\frac{1}{2}{\tilde{\mathit{\omega }}}^{\text{T}}\dot{M}(\mathit{q})\tilde{\mathit{\omega }}\end{array}$

which along the trajectories of the nominal system in equation (18) yields

$\begin{array}{c}\dot{V}(\tilde{\mathit{q}},\mathit{\xi },\tilde{\mathit{\omega }})={\tilde{\mathit{q}}}^T\tilde{\mathit{\omega }}-{\tilde{\mathit{q}}}^T\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\tilde{\mathit{q}}-{\alpha }^2{\mathit{\xi }}^T{K}_{vp}\mathit{\xi }\\ -{\tilde{\mathit{\omega }}}^T{K}_{vp}\tilde{\mathit{\omega }}-2\alpha \mathit{\gamma }{\mathit{\xi }}^T\tilde{\mathit{\omega }}\end{array}$

23

where the Property 2 has been used, as well as equations (12) and (13). The upper bounds of each term of equation (23) can be written as

$\begin{array}{c}{\tilde{\mathit{q}}}^T\tilde{\mathit{\omega }}\hspace{1em}\le \hspace{1em}\parallel \tilde{\mathit{q}}\parallel \parallel \tilde{\mathit{\omega }}\parallel ,\\ -\hspace{0.17em}{\tilde{\mathit{q}}}^T\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\hspace{0.17em}\tilde{\mathit{q}}\hspace{1em}\le \hspace{1em}-{\mathit{\lambda }}_{\Phi l}\parallel \tilde{\mathit{q}}{\parallel }^2,\\ -\hspace{0.17em}{\alpha }^2{\xi }^T{K}_{vp}\xi \hspace{1em}\le \hspace{1em}-{\alpha }^2(k_p+\mathit{\gamma })\parallel \xi {\parallel }^2,\\ -\hspace{0.17em}{\tilde{\mathit{\omega }}}^T{K}_{vp}\tilde{\mathit{\omega }}\hspace{1em}\le \hspace{1em}-(k_p+\mathit{\gamma })\parallel \tilde{\mathit{\omega }}{\parallel }^2,\\ -\hspace{0.17em}2\alpha \mathit{\gamma }{\xi }^T\tilde{\mathit{\omega }}\hspace{1em}\le \hspace{1em}2\alpha \mathit{\gamma }\parallel \xi \parallel \parallel \tilde{\mathit{\omega }}\parallel \end{array}$

where ${\mathit{\lambda }}_{\Phi l}$ is defined as

${\mathit{\lambda }}_{\Phi l}=\underset{\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}}\in {ℝ}^n}{{\displaystyle \text{inf}}}\left({\mathit{\lambda }}_{\text{min}}\left\{\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})\right\}\right)$

Thus, equation (23) can be upper bounded as

$\dot{V}(\tilde{\mathit{q}},\xi ,\tilde{\mathit{\omega }})\le -{\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]}^{\text{T}}Q\left[\begin{array}{c}\parallel \tilde{\mathit{q}}\parallel \\ \parallel \xi \parallel \\ \parallel \tilde{\mathit{\omega }}\parallel \end{array}\right]$

24

where

$Q=\left[\begin{array}{ccc}{\mathit{\lambda }}_{\Phi l}& 0& -\frac{1}{2}\\ 0& {\alpha }^2(k_p+\mathit{\gamma })& -\alpha \mathit{\gamma }\\ -\frac{1}{2}& -\alpha \mathit{\gamma }& k_p+\mathit{\gamma }\end{array}\right]$

To ensure that $\dot{V}(t,\mathit{z})$ is negative definite, matrix Q has to be positive definite, which is satisfied if

${\mathit{\lambda }}_{\Phi l}>\frac{k_p+\mathit{\gamma }}{4k_p\left(k_p+2\mathit{\gamma }\right)}$

25

From equation (24), $\dot{V}(t,\mathit{z})$ can also be bounded as

$\dot{V}(t,\mathit{z})\le -{\mathit{\lambda }}_3\parallel \mathit{z}{\parallel }^2$

26

where

${\mathit{\lambda }}_3={\mathit{\lambda }}_{\text{min}}\left\{Q\right\}$

27

Thus, it can be concluded the exponential stability of the origin of the nominal system by satisfying equations (22) and (25).

Stability of the overall closed-loop system

The temporal derivative of the Lyapunov candidate function along the trajectories of the overall closed-loop system in equation (17), after simplifying terms can be written as

$\begin{array}{l}\dot{V}(\tilde{q},\xi ,\tilde{\omega })={\tilde{q}}^{\text{T}}\tilde{\omega }-{\tilde{q}}^{\text{T}}\Phi (\tilde{q},\dot{\tilde{q}})\hspace{0.17em}\tilde{q}-{\alpha }^2{\xi }^{\text{T}}{K}_{vp}\xi \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\tilde{\omega }}^{\text{T}}{K}_{vp}\tilde{\omega }-2\alpha \gamma {\xi }^{\text{T}}\tilde{\omega }+({\tilde{\omega }}^{\text{T}}+\alpha {\xi }^{\text{T}})\hspace{0.17em}\mathit{d}(t,\mathit{z})\end{array}$

28

The last term of equation (28) can be upper bounded as

$({\tilde{\mathit{\omega }}}^{\text{T}}+\alpha {\mathit{\xi }}^{\text{T}})\hspace{0.17em}\mathit{d}\hspace{0.17em}(t,\mathit{z})\le (1+\alpha )\parallel \mathit{z}\parallel \parallel \mathit{d}(t,\mathit{z})\parallel$

29

where it has been taken into account that $\parallel \tilde{\mathit{\omega }}\parallel \le \parallel \mathit{z}\parallel$ , $\parallel \mathit{\xi }\parallel \le \parallel \mathit{z}\parallel$ . The perturbation can be upper bounded as (see equation (1U) in Appendix 1)

$\parallel \mathit{d}(t,\mathit{z})\parallel \le {\varsigma }_0+{\varsigma }_1\parallel \mathit{z}\parallel +{\varsigma }_2{\parallel \mathit{z}\parallel }^2$

30

where ${\varsigma }_0$ , ${\varsigma }_1$ , and ${\varsigma }_2$ are positive constants (see equations (1V) to (1X) in Appendix 1). By substituting equation (30) in equation (29), we have

$\begin{array}{l}\left({\tilde{\omega }}^{\text{T}}+\alpha {\xi }^{\text{T}}\right)\hspace{0.17em}\mathit{d}\hspace{0.17em}(t,\mathit{z})\le \left(1+\alpha \right){\varsigma }_0\parallel \mathit{z}\parallel +\left(1+\alpha \right){\varsigma }_1{\parallel \mathit{z}\parallel }^2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(1+\alpha \right){\varsigma }_2{\parallel \mathit{z}\parallel }^3\end{array}$

31

After substituting equations (26) and (31) as the upper bounds of equation (28) and making some arrangements, we obtain

$\begin{array}{l}\dot{V}\le (1+\alpha ){\varsigma }_0\parallel \mathit{z}\parallel -({\lambda }_3-(1+\alpha ){\varsigma }_1){\parallel \mathit{z}\parallel }^2\\ \hspace{1em}\hspace{1em}+(1+\alpha ){\varsigma }_2{\parallel \mathit{z}\parallel }^3\end{array}$

32

Now we define a ball ${\mathcal{B}}_r\subset {ℝ}^n$ , centered in the origin, with radio r, such that

${\mathcal{B}}_r:=\left\{\mathit{z}\in {ℝ}^n:\parallel \mathit{z}\parallel <r\right\}$

33

inside which there exist conditions for the temporal derivative of the Lyapunov candidate function $\dot{V}$ to be negative definite, as it will be addressed later. Equation (32) can now be written as

$\dot{V}\le (1+\alpha ){\varsigma }_0\parallel \mathit{z}\parallel -({\mathit{\lambda }}_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2)){\parallel \mathit{z}\parallel }^2$

34

Let introduce a positive constant $\epsilon <1$ , such that equation (34) now can be written as

$\begin{array}{l}\dot{V}\le (1+\alpha ){\varsigma }_0\parallel \mathit{z}\parallel \\ \hspace{1em}\hspace{1em}-({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2))(1-\epsilon ){\parallel \mathit{z}\parallel }^2\\ \hspace{1em}\hspace{1em}-({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2))\epsilon {\parallel \mathit{z}\parallel }^2\end{array}$

35

Investigating the positiveness of the sum of the first and the third terms of equation (35) we have that

$-({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2))\epsilon \parallel \mathit{z}{\parallel }^2+(1+\alpha ){\varsigma }_0\parallel \mathit{z}\parallel <0$

when

$\begin{array}{l}({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2))\epsilon \parallel \mathit{z}{\parallel }^2>(1+\alpha ){\varsigma }_0\parallel \mathit{z}\parallel ,\\ \parallel \mathit{z}\parallel \hspace{0.17em}>\hspace{0.17em}\frac{(1+\alpha ){\varsigma }_0}{({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2))\epsilon }\end{array}$

Thus, by denoting

$\mu :=\frac{\left(1+\alpha \right){\varsigma }_0}{\left({\lambda }_3-(1+\alpha )({\varsigma }_1+r{\varsigma }_2)\right)\epsilon }$

36

from equation (35) we can write

$\dot{V}\le -\left({\lambda }_3-\left(1+\alpha \right)\left({\varsigma }_1+r{\varsigma }_2\right)\right)\left(1-\epsilon \right)\parallel \mathit{z}{\parallel }^2,\hspace{0.17em}\forall \hspace{0.17em}\parallel \mathit{z}\parallel >\mu$

37

Notice that the coefficient of $\parallel \mathit{z}{\parallel }^2$ in equation (37) is negative, and the temporal derivative of the Lyapunov candidate function $\dot{V}$ is negative definite whenever

${\mathit{\lambda }}_3>\left(1+\alpha \right)\left({\varsigma }_1+r{\varsigma }_2\right)\left(1-\epsilon \right)$

38

is satisfied.

Now we are ready to establish the main stability result in the following.

Proposition 1

Let consider the dynamic model of a parallel robot in equation (7), together with the control law given by equations (9) to (11), where $\Phi (\tilde{\mathit{q}},\dot{\tilde{\mathit{q}}})$ is a diagonal matrix of variable gains, and $K_{vp}$ and $K_{vi}$ are diagonal matrices of constant gains defined in equations (12) to (13). Also, the parameter $k_p$ satisfies equation (22). In addition, the parameter μ in equation (36), where ${\mathit{\lambda }}_3$ was defined in equation (27), ${\varsigma }_0$ , ${\varsigma }_1$ , and ${\varsigma }_2$ are defined in equations (1V) to (1X) (see Appendix 1), satisfies

$\mu <\sqrt{\frac{{\mathit{\lambda }}_1}{{\mathit{\lambda }}_2}}\upsilon$

where ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ were defined in equation (21), and υ is the radius of a ball $\mathcal{D}:=\mathit{z}\in {ℝ}^{3n}\hspace{0.17em}:\hspace{0.17em}\parallel \mathit{z}\parallel <\upsilon$ , inside which the temporal derivative of the Lyapunov function in equation (37) is negative definite. Moreover, ${\mathit{\lambda }}_3$ satisfies equation (38). Then, the solutions $\mathit{z}(t)$ of the closed-loop system in equation (17) with initial state

$\mathit{z}(t_0)<\sqrt{\frac{{\mathit{\lambda }}_1}{{\mathit{\lambda }}_2}}\upsilon$

39

are uniformly ultimately bounded with ultimate bound given by

$\parallel \mathit{z}\parallel \le \sqrt{\frac{{\mathit{\lambda }}_2}{{\mathit{\lambda }}_1}}\mu ,\hspace{0.17em}\forall \hspace{0.17em}t>t_0+T$

with $T\ge 0$ . A Lyapunov function to demonstrate this is equation (19).

Proof

First, note in equation (20) that the Lyapunov candidate function V is upper and lower bounded by positive definite functions of $\parallel \mathit{z}\parallel$ . Moreover, the time derivative of V along the trajectories of the closed-loop system is upper bounded by a negative definite function of $\parallel \mathit{z}\parallel$ (see equation (37)), for all $\parallel \mathit{z}\parallel \ge \mu$ where μ is a positive parameter. By invoking the Theorem 4.18 in Khalil’s study ³² and recalling the described functions and parameters, the Proposition 1 is proven.

Remark 1

For any initial state $\mathit{z}(t_0)<\sqrt{{\mathit{\lambda }}_1/{\mathit{\lambda }}_2}\upsilon$ (inside the ball $\mathcal{D}$ ) in equation (39), the trajectories of the solution $\mathit{z}(t)$ of the system will stay below the uniform ultimate bound $\sqrt{{\mathit{\lambda }}_2/{\mathit{\lambda }}_1}\mu$ , after a finite time $t_0+T$ . This bound can be reduced by adjusting the bounds of the perturbation, which depends on the bounds of the desired trajectory and on the bounds of the variable gains (see equations (30) and equations (1V) to (1X) in Appendix 1).

Remark 2

Notice that the radius υ of the ball $\mathcal{D}$ is not necessarily equal to the radius r of the ball ${\mathcal{B}}_r$ in equation (33). A further analysis of this can be consulted in the study by Salas et al. ³³

Simulations

To test the effectiveness of the proposed control approach, we conducted numerical simulations of the proposed controller applied to the model of a planar parallel 3-revolute-revolute-revolute (RRR) joints robot. In the following subsection, the matrices of the dynamic model of the parallel 3-RRR robot are presented. Next, the desired trajectory for the end-effector of the robot is described. Finally, the simulations conducted are described.

Kinematic structure of the parallel 3-RRR robot

The kinematic structure of the robot is shown in Figure 5. This structure is designated as 3-RRR, which represents a robot with three legs or kinematic chains, each having three rotational joints, from which the joint to the fixed platform is actuated. ³³ The joints to the fixed platform are located in the dots labeled $A_1$ , $A_2$ , and $A_3$ . The mobile platform, which is considered the end-effector of the robot, is a triangular equilateral platform with $P=\left({\mathit{x}}_p,y_p\right)$ as its geometric center.

Figure 5.

Structure of the parallel robot 3-RRR.

The equation of the dynamic model of the robot used for simulations is equation (1). The nonzero elements of the matrix $M(\mathit{\rho })\in {ℝ}^{9\times 9}$ are ²⁹

$\begin{array}{c}{M}_{11}=M_{22}=M_{33}=s_1,\\ M_{44}=M_{55}=M_{66}=s_2,\\ M_{77}=M_{88}=M_{99}=m_p,\\ M_{14}=M_{41}=s_3\text{cos}({\mathit{q}}_1-{\beta }_1),\\ M_{25}=M_{52}=s_3\text{cos}({\mathit{q}}_2-{\beta }_2),\\ M_{36}=M_{63}=s_3\text{cos}({\mathit{q}}_3-{\beta }_3)\end{array}$

where $s_1=m_1{l}_{c1}^2+m_2{l}_1^2+I_1$ , $s_2=m_2{l}_{c2}^2$ , $s_3=m_2{l}_1{l}_{c2}$ .

The nonzero elements of the matrix $C\left(\mathit{\rho },\dot{\mathit{\rho }}\right)\in {ℝ}^{9\times 9}$ are

$\begin{array}{l}{C}_{14}=-s_4\text{sin}\left({\mathit{q}}_1-{\beta }_1\right){\dot{\mathit{q}}}_1,\\ C_{41}=s_4\text{sin}\left({\mathit{q}}_1-{\beta }_1\right){\dot{\beta }}_1,\\ C_{25}=-s_4\text{sin}\left({\mathit{q}}_2-{\beta }_2\right){\dot{\mathit{q}}}_2,\\ C_{52}=s_4\text{sin}\left({\mathit{q}}_2-{\beta }_2\right){\dot{\beta }}_2,\\ C_{36}=-s_4\text{sin}\left({\mathit{q}}_3-{\beta }_3\right){\dot{\mathit{q}}}_3,\\ C_{63}=s_4\text{sin}\left({\mathit{q}}_3-{\beta }_3\right){\dot{\beta }}_3\end{array}$

where $s_4=m_2{l}_1{l}_{c2}$ .

Since the robot workspace is limited to the horizontal plane, it can be considered that the robot is not affected by gravity forces. The nonzero elements of the matrix $D(\mathit{\rho })=\partial \mathit{\gamma }(\mathit{\rho })/\partial \mathit{\rho }\in {ℝ}^{(6)\times (9)}$ are

$\begin{array}{c}{D}_{11}=-l_1\text{sin}({\mathit{q}}_1),\\ D_{14}=-l_2\text{sin}({\beta }_1),\\ D_{19}=-r\text{sin}(\varphi +{\varphi }_1),\\ D_{21}=l_1\text{cos}({\mathit{q}}_1),\\ D_{24}=l_2\text{cos}({\beta }_1),\\ D_{29}=r\text{cos}(\varphi +{\varphi }_1),\\ D_{32}=-l_1\text{sin}(q_2),\\ D_{35}=-l_2\text{sin}({\beta }_2),\\ D_{39}=-r\text{sin}(\varphi +{\varphi }_2),\\ D_{42}=l_1\text{cos}({\mathit{q}}_2),\\ D_{45}=l_2\text{cos}({\beta }_2),\\ D_{49}=r\text{cos}(\varphi +{\varphi }_2),\\ D_{53}=-l_1\text{sin}({\mathit{q}}_3),\\ D_{56}=-l_2\text{sin}({\beta }_3),\\ D_{59}=-r\text{sin}(\varphi +{\varphi }_3),\\ D_{63}=l_1\text{cos}({\mathit{q}}_3),\\ D_{66}=l_2\text{cos}({\beta }_3),\\ D_{69}=r\text{cos}(\varphi +{\varphi }_3)\end{array}$

where ${\varphi }_1=\pi /6$ , ${\varphi }_2=5\pi /6$ , and ${\varphi }_1=3\pi /2$ [radians]. The vector of Lagrange multipliers $\mathit{\lambda }\in {ℝ}^6$ can be computed by the expression

$\begin{array}{l}\lambda ={\left[D(\rho )M^{-1}(\rho )D^{\text{T}}(\rho )\right]}^{-1}\\ \hspace{1em}•\hspace{0.17em}\left[D(\rho )M^{-1}(\rho )(C(\rho ,\dot{\rho })\dot{\rho }-\tau )-\dot{D}(\rho ,\dot{\rho })\dot{\rho }\right]\end{array}$

(see Garcia-Gamez et al.’s work ²⁹ ).

The parameters of the dynamic model are shown in Table 2. The inverse kinematic model ³⁴ is

${\mathit{q}}_i=2\hspace{0.17em}{{\displaystyle \text{tan}}}^{-1}\frac{-E_{1i}\pm \sqrt{E_{1i}^2+E_{2i}^2-E_{3i}^2}}{E_{3i}-E_{2i}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,3$

where

$E_{11}=-2l\left(y_p-\frac{L\text{sin}\varphi }{2}-\frac{\sqrt{3}L\text{cos}\varphi }{6}\right)$

40

$E_{21}=-2l\left({\mathit{x}}_p-\frac{L\text{cos}\varphi }{2}+\frac{\sqrt{3}L\text{sin}\varphi }{6}\right)$

41

$\begin{array}{l}{E}_{31}={\left({\mathit{x}}_p-\frac{L\text{cos}\varphi }{2}+\frac{\sqrt{3}L\text{sin}\varphi }{6}\right)}^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{1em}+{\left(y_p-\frac{L\text{sin}\varphi }{2}-\frac{\sqrt{3}L\text{cos}\varphi }{6}\right)}^2\end{array}$

42

$E_{12}=-2l\left(y_p+\frac{L\text{sin}\varphi }{2}-\frac{\sqrt{3}L\text{cos}\varphi }{6}\right)$

43

$E_{22}=-2l\left({\mathit{x}}_p+\frac{L\text{cos}\varphi }{2}+\frac{\sqrt{3}L\text{sin}\varphi }{6}-R\right)$

44

$E_{32}={\left({\mathit{x}}_p+\frac{L\text{cos}\varphi }{2}+\frac{\sqrt{3}L\text{sin}\varphi }{6}-R\right)}^2+{\left(y_p+\frac{L\text{sin}\varphi }{2}-\frac{\sqrt{3}L\text{cos}\varphi }{6}\right)}^2$

45

$E_{13}=-2l\left(y_p+\frac{\sqrt{3}L\text{cos}\varphi }{3}-\frac{\sqrt{3}}{2}R\right)$

46

$E_{23}=-2l\left({\mathit{x}}_p-\frac{\sqrt{3}L\text{sin}\varphi }{3}-\frac{1}{2}R\right)$

47

$E_{33}={\left({\mathit{x}}_p-\frac{\sqrt{3}L\text{sin}\varphi }{3}-\frac{1}{2}R\right)}^2+\hspace{0.17em}{\left(y_p+\frac{\sqrt{3}L\text{sin}\varphi }{3}-\frac{\sqrt{3}}{2}R\right)}^2$

48

with $L=0.5\hspace{0.17em}\text{m}$ as the length of the side of the mobile platform, and $R=0.62\hspace{0.17em}\text{m}$ as the length of the fixed platform, from $A_1$ to $A_2$ (see Figure 5).

Table 2.

Dynamic parameters of the model.

Parameter	Value	Units
$m_1$	2.1992	kg
$m_2$	2.0485	kg
$m_p$	5.8579	kg
$I_1$	0.0264	kg $\cdot {\text{m}}^2$
$I_2$	0.0228	kg $\cdot {\text{m}}^2$
$I_p$	0.2504	kg $\cdot {\text{m}}^2$
l	0.265	m
$l_{c1}$	0.1365	m
$l_{c2}$	0.1463	m
r	0.2887	m

Desired trajectory

The desired trajectory for the end-effector of the manipulator is an elliptical trajectory in operational space is given by

$\begin{array}{c}{\mathit{x}}_d={\mathit{x}}_h+U_{\mathit{x}}\text{cos}\hspace{0.17em}(wt)\\ y_d=y_h+U_y\text{sin}\hspace{0.17em}(wt)\end{array}$

where ${\mathit{x}}_h$ and $y_h$ are the initial position coordinates of the end-effector. The desired trajectory in joint space is obtained by computing the inverse kinematic model given in equations (40) to (48). The values of the initial position used in simulations are: ${\mathit{x}}_h=0.31\hspace{0.17em}\text{m}$ and $y_h=0.1789\hspace{0.17em}\text{m}$ . The initial value of the orientation of the end-effector is ${\varphi }_h=1.4811\hspace{0.17em}\text{rad}$ , which is a constant value for the whole trajectory. This initial position yields the initial joint position vectors ${\mathit{q}}_h=[0.4033\hspace{0.17em}\hspace{0.17em}2.4984\hspace{0.17em}\hspace{0.17em}-{1.6903]}^{\text{T}}$ and ${\beta }_h=[-0.7835\hspace{0.17em}\hspace{0.17em}1.3101\hspace{0.17em}\hspace{0.17em}-{2.8781]}^{\text{T}}$ rad. The amplitudes used in simulations are $U_{\mathit{x}}=0.08\hspace{0.17em}\text{m}$ and $U_y=0.06\hspace{0.17em}\text{m}$ . The trajectory must be completed in a period of time of 8 s. The coordinates of the joints $A_1$ , $A_2$ , and $A_3$ (see Figure 5) are $A_1=[0\hspace{0.17em}\hspace{0.17em}0]$ , $A_2=[0.62\hspace{0.17em}\hspace{0.17em}0]$ , and $A_3=[0.31\hspace{0.17em}\hspace{0.17em}0.5369]\hspace{0.17em}\text{m}$ .

An important issue regarding the motion of parallel manipulators is the avoidance of singular configurations. These configurations exist in the workspace of the manipulators and are represented by the mathematical singularity of the jacobian matrices of the differential kinematic model. It is well known that passing or approaching a singular configuration in the trajectory of the robot can result in high joint velocities and important uprising of the mechanical load, including an eventual loss of the control of the mechanism. ³⁵ There exists, in the literature of motion of parallel robot, several approaches for singularity avoidance, from trajectory planning to redundancy schemes. ^36,37 In this work, we address this problem by designing a desired trajectory that is free of singular configurations. With that objective in mind, we compute the reciprocal of the condition number ${\Gamma }^{-1}(\Lambda )={\mathit{\sigma }}_{\mathit{\lambda }}/{\mathit{\sigma }}_{\Lambda }$ , as is proposed by Alba, ³⁷ where ${\mathit{\sigma }}_{\mathit{\lambda }}$ and ${\mathit{\sigma }}_{\Lambda }$ denote the minimum and maximum singular values of a jacobian matrix Λ, respectively. Notice that a ${\Gamma }^{-1}(\Lambda )=0$ corresponds to a singular configuration. The jacobian matrix computed is the matrix $J_{\mathit{x}}\left(\mathit{\rho }\right)$ of the differential kinematic model $\dot{\mathit{x}}=J_{\mathit{x}}(\mathit{\rho })\dot{\mathit{q}}$ . As a result of the computing of ${\Gamma }^{-1}(J_{\mathit{x}}(\mathit{\rho }))$ along the desired trajectory, the minimum computed value is 0.1117, which indicates that it never reaches a zero value.

Simulations

A practical issue concerning the actuators of a robot is to avoid saturating the actuators with torque control references that exceed the rated values, particularly at the beginning of the test. In this work, to avoid actuator saturation in a similar way to real-time experiments, we established an upper bound for the demanded torques from each actuator. This upper bound prevent the controllers to produce a torque control signal that a real actuator cannot accomplish in real-time experiments. The upper bound for the torque control signals was set to 4 Nm.

The selected gains of the inner loop PI control loop are shown in Table 3. These gains were selected by a trial and error process to avoid exceeding the upper bound of the actuators, particularly at the beginning of the test. For comparison purposes, simulations of two more controllers: a classic PID controller and a FPD controller applied to the 3-RRR parallel robot model were performed. The gains of the PID controller were chosen by means of a trial and error procedure, to obtain the best possible performance from the robot. Moreover, we were careful in avoiding to exceed the upper bound of the nominal torque of the actuators at the beginning of the test. These gains are shown in Table 4.

Table 3.

Gains of the PI control loop.

Parameter	Value	Units
$K_{vp}$	$\text{diag}\left\{5,10,8\right\}$	$\text{Nm}\cdot \text{s/rad}$
$K_{vi}$	$\text{diag}\left\{50,20,70\right\}$	$\text{Nm/rad}$

PI: proportional integral.

Table 4.

Gains of the PID controller.

Parameter	Value	Units
K p	diag {9, 50, 450}	Nm/rad
K i	diag {10, 10, 10}	Nm/rad·s
K d	diag {40, 70, 10}	Nm·s/rad

PID: proportional integral derivative.

The FPD controller is the proposed by Salas et al., ³³ which is inspired on the work by Llama et al. ³⁸ Both gains P and D are tuned by a fuzzy inference engine. Three input fuzzy membership functions and three output fuzzy membership functions are used. The partitions of these functions are shown in Table 5. Similar to the work reported by Salas et al., ³³ the inference mechanism is supported on three fuzzy inference rules

If ${\mu }_S(|{\tilde{\mathit{q}}}_i|)$ then ${\mu }_B(y)$

If ${\mu }_M(|{\tilde{\mathit{q}}}_i|)$ then ${\mu }_M(y)$

If ${\mu }_B(|{\tilde{\mathit{q}}}_i|)$ then ${\mu }_S(y)$

where ${\mu }_S(|{\tilde{\mathit{q}}}_i|)$ , ${\mu }_M(|{\tilde{\mathit{q}}}_i|)$ , and ${\mu }_M(|{\tilde{\mathit{q}}}_i|)$ represent the membership grades of the input variable ${\tilde{\mathit{q}}}_i$ to the fuzzy sets represented by linguistic variables Small (S), Medium (M), and Big (B). On the other hand, ${\mu }_B(y)$ , ${\mu }_M(y)$ , and ${\mu }_S(y)$ represent the fuzzy membership grades to obtain the crisp variable from the output membership functions.

Table 5.

Partitions of the fuzzy membership functions of the FPD controller.

Joint	$k_{1p}$	$k_{2p}$	$k_{3p}$	$k_{1d}$	$k_{2d}$	$k_{3d}$
1	9.5	90	200	3.7	16.1	60
2	9.5	90	200	3.7	16.1	60
3	9.5	90	200	0.61	2.83	5.29
Units	[Nm/rad]			[Nm·s/rad]

FPD: fuzzy-tuned proportional derivative.

The defuzzification is carried out by a simplified version of the Center-of-Gravity method, ³⁸ which can be written as

$y=\frac{{\mathit{k}}^{\text{T}}\mu (|{\tilde{\mathit{q}}}_i|)}{\parallel \mu (|{\tilde{\mathit{q}}}_i|)\parallel }$

where $\mathit{k}={\left[k_1\hspace{0.17em}{k}_2\hspace{0.17em}{k}_3\right]}^{\text{T}}$ is the output partition.

Simulation results are shown in Figures 6 to 12. In Figures 6 to 8, the errors in joints 1 to 3, respectively, obtained by the three controllers, are depicted. The largest tracking error peaks obtained by the PID controller were 0.2, 0.18, and 0.0055 rad for the joints 1, 2, and 3, respectively, while the largest tracking error peaks obtained by the SOF-P-PI controller were 0.01, 0.01, and 0.005 rad, for the joints 1, 2, and 3, respectively. The largest tracking error peaks obtained by the FPD were 0.06, 0.02, and 0.0053 rad, for the joints 1, 2, and 3, respectively. It can be noticed in Figures 6 to 7 that both errors produced by the PID controller in joints 1 and 2 have much larger overshooting and settling times that the errors produced by the SOF-P-PI controller. On the other hand, the overshooting and settling times produced by the FPD controller are only a little larger than in the SOF-P-PI controller. In addition, it can be noticed in Figure 8 that although the overshooting response in the three controllers for the joint 3 are similar, the steady-state error of the SOF-P-PI controller is smaller than the steady-state error of the PID and the FPD controllers. In Figures 9 to 11 the demanded torques by the SOF-P-PI, the PID, and the FPD controllers, respectively, are shown. Observe in these figures that, in the entire simulation time, the torques are below the predefined value for maximum torques. In Figure 12, the temporal evolution of the variable P gains in the outer loop are shown. This temporal evolution exhibits the increase of the gain values according to the adopted criteria in the SOF algorithm, to obtain small tracking errors. Also, in Figures 13 and 14, the temporal evolution of the P and D gains, respectively, of the FPD controller are shown.

Figure 6.

Position errors in joint 1. SOF: self-organizing fuzzy; FPD: fuzzy-tuned proportional derivative; PID: proportional integral derivative; P-PI: proportional–proportional integral.

Figure 7.

Position errors in joint 2. SOF: self-organizing fuzzy; FPD: fuzzy-tuned proportional derivative; PID: proportional integral derivative; P-PI: proportional–proportional integral.

Figure 8.

Position errors in joint 3. SOF: self-organizing fuzzy; FPD: fuzzy-tuned proportional derivative; PID: proportional integral derivative; P-PI: proportional–proportional integral.

Figure 9.

Commanded torques by SOF-P-PI. SOF: self-organizing fuzzy; P-PI: proportional–proportional integral.

Figure 10.

Commanded torques by PID. PID: proportional integral derivative.

Figure 11.

Commanded torques by FPD. FPD: fuzzy-tuned proportional derivative.

Figure 12.

Evolution of P gains by SOF algorithm. P: proportional; SOF: self-organizing fuzzy.

Figure 13.

Evolution of P gains by FPD. P: proportional; FPD: fuzzy-tuned proportional derivative.

Figure 14.

Evolution of D gains by FPD. D: derivative; FPD: fuzzy-tuned proportional derivative.

A useful quantitative measurement of the performance of a controller can be done by computing the ${\mathcal{L}}_2$ norm of the position or tracking errors. The ${\mathcal{L}}_2$ norm can be obtained with the formula

${\parallel \tilde{\mathit{q}}\parallel }_{{\mathcal{L}}_2}=\sqrt{\frac{1}{T-t_0}{\displaystyle {\int }_{t_0}^T}\hspace{0.17em}{\tilde{\mathit{q}}}^{\text{T}}\tilde{\mathit{q}}\text{d}t}$

Here, T is the computed time period and $t_0$ is the initial computing time. The total computing time considered must be the period of steady state behavior, this is, it must start after the transient behavior has passed. The result of computing the ${\mathcal{L}}_2$ norm of the position errors from the three controllers, SOF-P-PI, PID, and FPD is shown in Table 6. The period of steady state behavior is from 2.5 s to 8 s.

Table 6.

${\mathcal{L}}_2$ norm of $\tilde{\mathit{q}}$ from the three controllers.

SOF-P-PI	PID	FPD
0.0022	0.1107	0.0245

SOF: self-organizing fuzzy; PID: proportional integral derivative; FPD: fuzzy-tuned proportional derivative; P-PI: proportional–proportional integral.

Conclusion

In this article, we have presented a novel two-loop P-PI tracking controller for parallel manipulators, in which the gains of the outer loop are online tuned by an original SOF algorithm. This controller is a not model-based scheme which improves the robustness in the presence of perturbations and parameter uncertainty of the constant gains two-loop approach for tracking control, by implementing a variable gains loop.

The stability analysis of the closed-loop system has been developed, resulting in the uniform ultimate boundedness of the solutions of the closed-loop system. It is worth to mention that this result can be generalized not only to P-PI tracking controllers with variable gains tuned by an SOF algorithm but also to other variable gains P-PI tracking control schemes for parallel robots, only providing that the gains produced by the tuning algorithm satisfy the boundedness criterion. This is an original result not published elsewhere.

The usefulness of the proposed controller was validated by performing numerical simulations of the application of the controller to a parallel manipulator. In addition, simulations of a PID controller and a FPD controller applied to the parallel manipulator have been conducted for comparison purposes. The simulation results of the three controllers were contrasted, from which a significant quantitative measurement result was obtained that shows the better results of the SOF-P-PI controller.