**About the Editors**

#### **Marco Carricato**

Prof. Marco Carricato received an M.Sc. degree (with honors) in Mechanical Engineering in 1998, and a Ph.D. degree in Mechanics of Machines in 2002. He has been working with the University of Bologna since 2004. He is currently a Full Professor in the Department of Industrial Engineering, where he is the head of the IRMA L@B (Industrial Robotics, Mechatronics & Automation Lab @ Bologna). His research interests include robotic systems, servo-actuated automatic machinery and the theory of mechanisms, with a particular emphasis on parallel manipulators (displacement analysis, kinematics, dynamics, synthesis, gravity compensation, cable drives), mobile collaborative robotics, efficiency and optimization of servomechanisms, and the theory of 'screws'. His work in the aforementioned areas has been the subject of a number of scientific publications in international conferences and journals.

#### **Edoardo Idà**

Dr. Edoardo Idà received B.Sc. and M.Sc. degrees (with honors) in Mechanical Engineering at the University of Bologna, and a Ph.D. degree in Mechanics and Advanced Engineering Science at the same university in 2021. He is currently a Junior Assistant Professor in the Department of Industrial Engineering, where he is also the head of operations at IRMA L@B (Industrial Robotics, Mechatronics & Automation Lab @ Bologna). His research focuses on cable-driven robotic systems, continuum robots, automation, and mechanism design.

**Xiangquan Li 1,2 and Zhengguang Xu 1,\***


**Abstract:** This work addresses a pattern-moving-based partial form dynamic linearization model free adaptive control (P-PFDL-MFAC) scheme and illustrates the bounded convergence of its tracking error for a class of unknown nonaffine nonlinear discrete-time systems. The concept of pattern moving is to take the pattern class of the system output condition as a dynamic operation variable, and the control purpose is to ensure that the system outputs belong to a certain pattern class or some desired pattern classes. The P-PFDL-MFAC scheme mainly includes a modified tracking control law, a deviation estimation algorithm and a pseudo-gradient (PG) vector estimation algorithm. The classification-metric deviation is considered as an external disturbance, which is caused by the process of establishing the pattern-moving-based system dynamics description, and an improved cost function is proposed from the perspective of a two-player zero-sum game (TP-ZSG). The bounded convergence of the tracking error is rigorously proven by the contraction mapping principle, and the validity of the theoretical results is verified by simulation examples.

**Keywords:** pattern moving; partial form dynamic linearization (PFDL); nonlinear system; two-player zero-sum game; model free adaptive control (MFAC)

#### **1. Introduction**

In the process of industrial production, there is a range of complex equipment, such as sintering machines, rotary kilns, blast furnaces, and so on. Due to the increase in complexity, such as nonlinearity, high order, large delay, time-varying, and parameter perturbation, it is very difficult to establish an accurate mathematical model [1]. To a certain extent, this kind of production system is mainly governed by the law of statistical moving rather than the existing Newton's law of mechanics. A group of the same or similar system working conditions can produce the corresponding products with the same or similar quality index parameters [2].

A feasible method of system modeling and control is the pattern recognition technology for these considered systems [3], and most researchers' practice is to design the corresponding model and controller according to the different pattern classes of the system working condition [4,5]. Different from the previous multi-controller model design method based on pattern classes, a novel pattern-moving-based system dynamics description method was proposed in [6]. Its basic idea is to take the pattern class as a moving variable, and this variable is mapped to a computable space by class centers [7], interval numbers [8], and cells [9] due to its lack of arithmetic operation attribute. One advantage of the system dynamics description method introduced in [6] is that it is robust to system parameter disturbance and measurement noise. Regarding robust control, a well-known method is sliding mode control [10–12], which has a good ability to deal with external disturbances and system uncertainties. In recent years, a series of important research achievements have been made in sliding mode control, and many improved methods have been proposed,

**Citation:** Li, X.; Xu, Z. Pattern-Moving-Based Partial Form Dynamic Linearization Model Free Adaptive Control for a Class of Nonlinear Systems. *Actuators* **2021**, *10*, 223. https://doi.org/10.3390/act10090223

Academic Editor: Marco Carricato and Edoardo Idà

Received: 4 July 2021 Accepted: 2 September 2021 Published: 5 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

such as global sliding mode control [13] and terminal sliding mode control [14]. Different from the methods proposed in [10–14], the pattern-moving-based system dynamics description method is able to eliminate the system disturbance in the process of pattern classification, as long as the influence of the disturbance on the output does not change the pattern class to which the output belongs. In the case of various metric methods of pattern class, the linear autoregressive model with exogenous input (ARX) or interval ARX (IARX) model has been established, and the minimum-variance-based controller [6], optimal controller [15], predictive controller [16], and state-feedback-based [7] controller have been designed. However, it is well known that it is not easy to identify the system model order and parameters. In addition, even if a pattern-moving-based mathematical prediction model such as ARX or IARX is proposed, it is always an approximation of the real plant, and the unmodeled dynamics of the system are inevitable. Therefore, it is of significance to propose a pattern-moving-based data-driven control (DDC) method and design a controller whose parameters are adjusted by adopting the online input/output (I/O) data and the offline historical data simultaneously.

The data-driven controller is designed directly depending on the offline or/and online I/O data, instead of the explicit mathematical model of the controlled plant [17]. Generally, DDC can be almost cataloged into the following classes according to the different ways in which the data are used: (1) adaptive dynamic programming [18] and iterative learning control [19] based on offline and online data; (2) iterative feedback tuning [20] and virtual reference feedback tuning [21] based on offline data; (3) traditional MFAC [17,22–24] based on online data. The traditional MFAC method does not use the state space model but puts forward new concepts such as pseudo-gradient (PG) vector or pseudo-partial derivative (PPD) to capture the dynamic characteristics of the controlled plant, and it designs the controller through the dynamic linearization data model of the controlled plant at each operating point. Thus far, three equivalent dynamic linearization data models have been proposed, i.e., PFDL, compact-form dynamic linearization (CFDL), and the fullform dynamic linearization (FFDL) data model. By setting input correlation and output correlation components with different memory lengths, the three kinds of data models are different equivalent descriptions of system evolution, and they have different dynamic description capabilities for the controlled plant. Recently, due to many advantages of the MFAC method, such as the fact that establishing a controller merely depends on the measurement I/O data, the monotonic convergence of tracking error, and the boundedinput bounded-output stability of the closed-loop system, it has achieved many application results in many fields, and a few examples are as follows: the MFAC-based fault-tolerant control [25]; sensorless brushless direct current motor based on MFAC [26]; multi-agent systems tracking control [27]; MFAC-based sliding mode control [28]; chemical process based on MFAC [29], etc.

However, although the traditional MFAC algorithms have good control qualities for single-input single-output (SISO), multiple-input single-output, and multiple-input multiple-output time-varying structures and parameters in nonlinear discrete-time systems, there are few reports on MFAC for single-input multiple-output (SIMO) nonlinear systems or systems where the desired exact value of the output target cannot be determined exactly. In view of this kind of nonlinear system, a P-PFDL-MFAC method is proposed in this work and it considers that the difference in the output between next time and the current time is related to the differences in inputs in a time window between the current time and a specific previous time. The length of the time window corresponds to the number of PG vector elements, which is also called the pseudo-order of the equivalent PFDL data model. This is the most significant difference between the method proposed here and the pattern-moving-based CFDL-MFAC (P-CFDL-MFAC) scheme in [30], which considered that the output difference between next time and the current time is only related to the input difference between the current time and the previous time. The control purpose of this kind of system is to make the system outputs belong to one or some specific pattern classes. The first contribution of this work is to combine the pattern-moving-based system

dynamics description with the traditional PFDL-MFAC method, and to design a control law algorithm based on two-player zero-sum game and saddle point theory [31,32] under the condition of classification-metric deviation. Another major contribution is that the bounded convergence of the tracking error dynamics of the closed-loop control system is rigorously proven by using the contraction mapping principle.

The remainder of this work is organized as follows. Section 2 introduces the preliminary of the work. Section 3 presents the problem formulation and designs a patternmoving-based PFDL-MFAC scheme. The bounded convergence of the closed-loop system's tracking error is proven in Section 4. Section 5 presents two simulation examples to demonstrate the correctness and efficiency of the proposed algorithms. A conclusion is given in Section 6.

Notation: R denotes the real number domain; Z<sup>+</sup> denotes the positive integer domain; <sup>R</sup>*<sup>n</sup>* is the real *<sup>n</sup>*-dimensional space; [·] *<sup>T</sup>* is the transpose of [·]; -· is the Euclidean norm, and -·*v* is the consistent matrix norm.

#### **2. Preliminary**

Consider a class of SIMO nonaffine nonlinear discrete-time systems with unknown structure, order and parameters.

$$\begin{cases} y\_1(k+1) = f\_1(y\_1(k), \dots, y\_1(k-n\_1), \mu(k), \dots, \mu(k-m\_1)) + d\_1(k), \\ y\_2(k+1) = f\_2(y\_2(k), \dots, y\_2(k-n\_2), \mu(k), \dots, \mu(k-m\_2)) + d\_2(k), \\ \vdots \\ y\_q(k+1) = f\_q(y\_n(k), \dots, y\_n(k-n\_q), \mu(k), \dots, \mu(k-m\_q)) + d\_q(k), \end{cases} \tag{1}$$

where *<sup>q</sup>* <sup>&</sup>gt; 1; *yi*(*k*) denotes the output of *fi*(·) and it satisfies *yi*(*k*) <sup>∈</sup> <sup>R</sup>; *<sup>u</sup>*(*k*) is the whole system input and it satisfies *<sup>u</sup>*(*k*) <sup>∈</sup> <sup>R</sup>; *mi* , *ni* represent the unknown input and output orders, respectively, and they satisfy that *mi* <sup>∈</sup> <sup>Z</sup><sup>+</sup> , *ni* <sup>∈</sup> <sup>Z</sup>+; *di*(*k*) is the weak output measurement noise; *fi*(·) denotes an unknown nonlinear discrete-time function; *i* ∈ {1, ··· , *q*}.

**Assumption 1.** *The input of this kind of system (1) is bounded, i.e., a constant M*<sup>1</sup> *exists and satisfies that* |*u*(*k*)| ≤ *M*1*.*

A pattern-moving-based system dynamics description [6–9,30] that corresponds to system (1) is proposed in the following steps.


(3) Establishing the pattern-moving-based system dynamics equations. The inputs {*u*(*k*)}, implicit metric values { ¯*dx*(*k*)}, and class center explicit metric values {*s*(*k*)} are employed to construct the following dynamics equations.

$$
\bar{d}\mathbf{x}(k+1) = f(\bar{d}\mathbf{x}(k), \dots, \bar{d}\mathbf{x}(k-n), u(k), \dots, u(k-m)), \tag{2}
$$

$$s(k+1) = D(M(\bar{d}\mathbf{x}(k+1))) = \begin{cases} s\_{1\prime} \, \bar{d}\mathbf{x}(k+1) \in [s\_1 - r\_1, s\_1 + r\_1], \\ s\_{2\prime} \, \bar{d}\mathbf{x}(k+1) \in (s\_2 - r\_2, s\_2 + r\_2], \\ \vdots \\ s\_{N\prime} \, \bar{d}\mathbf{x}(k+1) \in (s\_N - r\_N, s\_N + r\_N], \end{cases} \tag{3}$$

where *f*(·) is an unknown SISO nonlinear discrete-time system function; *m*, *n* denote the input and output orders of system (2), respectively.

By choosing a reasonable classification method, such as a modified quantized control classification [34], it can be obtained that *Ci* = *si* + *ri* = *si*+<sup>1</sup> − *ri*+1, which is named the class threshold. It exits a classification-metric deviation *e*(*k* + 1) between the ¯*dx*(*k* + 1) and *<sup>s</sup>*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>), and <sup>|</sup>*e*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)<sup>|</sup> <sup>=</sup> <sup>|</sup>*s*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> ¯*dx*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)| ≤ *ri*, while *<sup>s</sup>*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) = *si*. Let *rmax* = *maxi*∈[1,*N*]{*ri*}, then |*e*(*k*)| ≤ *rmax*.

**Remark 1.** *As mentioned in the Introduction, the description of system dynamics based on pattern moving was first proposed in [6], and further studied in [7–9,30]. The basic idea is to treat the pattern class as a moving variable. Since this variable does not have the attribute of arithmetic operation, it is necessary to measure it into a computable space, and then construct the corresponding dynamic equation in this space. Obviously, the SISO nonlinear system or linear time-varying system can also be treated by the dynamic description method proposed in this section, but the feature extraction* (*T*(·)) *process is not required.*

**Remark 2.** *The ultimate goal of classifying and measuring the first principal component information is to obtain a SISO system dynamics description in a computable space. From the perspective of pattern recognition technology, when the contribution rate of the first principal component obtained after feature extraction is more than* 85%*, it is considered that the first principal component information does not lose the original information or it loses very little. If the contribution rate of the first principal component information does not reach* 85%*, more principal component information should be considered. Then, after classification and class center explicit metric, the metric result of each pattern class variable is a vector. A pattern-moving-based SIMO system dynamics description is to be constructed in a computable space, but the output dimension may be less than that of the original system. For the pattern-moving-based SIMO system, its control method remains to be studied in the future. In this work, we only consider the case in which the contribution rate of the first principal component information is greater than* 85%*.*

#### **3. Problem Formulation and Control Scheme**

#### *3.1. Problem Formulation*

Through the above system dynamics description method, the model free adaptive tracking control problem of system (1) is transformed into the corresponding control problem of system (2) and (3). In order to carry out our next analysis, the following assumptions and lemma are proposed first.

**Assumption 2.** *The partial derivatives of nonlinear system function f*(·) *with respect to all variables of the system (2) exist and are continuous.*

**Assumption 3.** *The system (2) satisfies the generalized Lipschitz condition, i.e.,*

$$|d\bar{\mathfrak{a}}(k\_1+1) - d\bar{\mathfrak{a}}(k\_2+1)| \le b || \mathcal{U}\_l(k\_1) - \mathcal{U}\_l(k\_2)||\_{\prime}$$

*where Ul*(*k*)=[*u*(*k*), ··· , *<sup>u</sup>*(*<sup>k</sup>* <sup>−</sup> *<sup>l</sup>* <sup>+</sup> <sup>1</sup>)]*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*<sup>l</sup> , l denotes the input pseudo-order, which satisfies l* > 1*, and b is a positive constant.*

**Lemma 1** ([22,23])**.** *For the considered system (2) satisfying Assumptions 2 and 3, there must exist a time-varying parameter vector ϕ<sup>f</sup>* ,*l*(*k*) *which is called a pseudo-gradient (PG) vector. If* -Δ*Ul*(*k*)-= 0*, the system (2) can be described as the following PFDL data model.*

$$
\Delta \bar{d} \mathbf{x}(k+1) = \varphi\_{f,l}^T(k) \Delta l I\_l(k), \tag{4}
$$

*where ϕ<sup>f</sup>* ,*l*(*k*)-<sup>≤</sup> *b;* <sup>Δ</sup> ¯*dx*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) = ¯*dx*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> ¯*dx*(*k*)*; <sup>ϕ</sup><sup>f</sup>* ,*l*(*k*)=[*ϕ*1(*k*), ··· , *<sup>ϕ</sup>l*(*k*)]*T;* Δ*Ul*(*k*) = *Ul*(*k*) − *Ul*(*k* − 1)*.*

Because the implicit metric values { ¯*dx*(*k*)} are not available, the traditional MFAC methods can not be directly used in such systems. Therefore, this work will focus on the design of a new control scheme that merely depends on the obtained data {*s*(*k*)}, {*u*(*k*)} and the performance analysis of the closed-loop control system.

#### *3.2. The P-PFDL-MFAC Scheme*

It can be seen from the system dynamics Equations (2) and (3) that there is a classificationmetric deviation *e*(*k* + 1) between the initial predicted output ¯*dx*(*k* + 1) and the final output *s*(*k* + 1) of the system, and this deviation *e*(*k* + 1) is always considered as a bounded external disturbance [12] in this work. Based on the saddle point theory of TP-ZSG proposed in [30–32], an improved cost function is designed in order to obtain a deviation estimation algorithm and an adaptive tracking control law, which aims to find an equilibrium point between the classification-metric deviation difference and the input difference. The basic idea is that even under large deviation fluctuation, a small input variation value can be found to optimize the loss function.

$$\begin{split} f(\Delta u(k), \Delta e(k+1)) &= \left| s^\*(k+1) - s(k+1) \right|^2 + \lambda \left| u(k) - u(k-1) \right|^2 \\ &- \gamma^2 \left| e(k+1) - e(k) \right|^2, \end{split} \tag{5}$$

where Δ*e*(*k* + 1) = *e*(*k* + 1) − *e*(*k*); *λ* is utilized to limit the variation in the control input difference, which satisfies *λ* > 0; *s*∗(*k* + 1) denotes the desired class center explicit metric value at time instant *k* + 1; *γ* is employed to limit the difference change in classificationmetric deviation, which satisfies *γ* > 1; Δ*u*(*k*) = *u*(*k*) − *u*(*k* − 1).

By solving the following equations

$$\frac{\partial f(\Delta u(k), \Delta e(k+1))}{\partial \Delta u(k)} = 0,$$

and

$$\frac{\partial f(\Delta u(k), \Delta e(k+1))}{\partial \Delta e(k+1)} = 0,$$

one has the optimal results, such as

$$
\Delta \epsilon (k+1) = \frac{1}{1 - \gamma^2} \left( s^\*(k+1) - s(k) - \boldsymbol{\varrho}\_{f,l}^T(k) \Delta l I\_l(k) \right), \tag{6}
$$

and

$$\Delta u(k) = \frac{\varrho\_1(k)\varrho\_1(s^\*(k+1) - s(k) - \Delta e(k+1))}{\lambda + |\varrho\_1(k)|^2} - \frac{\varrho\_1(k)\sum\_{i=2}^l \varrho\_i \varrho\_i(k)\Delta u(k-i+1)}{\lambda + |\varrho\_1(k)|^2}, \tag{7}$$

where *γ* > 1; *λ* > 0; *ρ<sup>i</sup>* is a step-size, which satisfies *ρ<sup>i</sup>* ∈ (0, 1] and makes the control algorithm more general; *i* ∈ {1, ··· , *l*}.

In order to estimate the PG vector, the following objective function is designed.

$$J(\boldsymbol{\varrho}\_{f,l}(k)) = \left| s(k) - s(k-1) - \boldsymbol{\varrho}\_{f,l}^T(k) \Delta l I\_l(k-1) \right|^2 + \mu \| \boldsymbol{\varrho}\_{f,l}(k) - \boldsymbol{\varrho}\_{f,l}(k-1) \|^2,\tag{8}$$

where *μ* is a weight factor and it satisfies *μ* > 0.

By letting

$$\frac{\partial J(\varphi\_{f,l}(k))}{\partial \varphi\_{f,l}(k)} = 0,$$

one can obtain the estimation algorithm of the PG vector as follows:

$$\varphi\_{f,l}(k) = \varphi\_{f,l}(k-1) + \frac{\eta \Delta l I\_l(k-1) \left(\Delta \mathbf{s}(k) - \varphi\_{f,l}^T(k-1)\Delta l I\_l(k-1)\right)}{\mu + \|\Delta l I\_l(k-1)\|^2},\tag{9}$$

where Δ*s*(*k*) = *s*(*k*) − *s*(*k* − 1); *η* is a step-size that satisfies *η* ∈ (0, 2] and makes the estimation algorithm more general; *μ* > 0.

Combining the above algorithms (6), (7), and (9), and proposing a reset algorithm of the PG estimation vector and a limitation mechanism of classification-metric deviation, the P-PFDL-MFAC scheme can be obtained.

$$\hat{\varphi}\_{f,l}(k) = \hat{\varphi}\_{f,l}(k-1) + \frac{\eta \Delta l I\_l(k-1) \left(\Delta s(k) - \hat{\varphi}\_{f,l}^T(k-1)\Delta l I\_l(k-1)\right)}{\mu + ||\Delta l I\_l(k-1)||^2},\tag{10}$$

$$\mathfrak{e}(k+1) = \mathfrak{e}(k) + \frac{1}{1 - \gamma^2} \left( \mathbf{s}^\*(k+1) - \mathbf{s}(k) - \mathfrak{q}\_{f,l}^T(k)\Delta l I\_l(k) \right),\tag{11}$$

$$\begin{split} u(k) &= u(k-1) - \frac{\phi\_1(k)\sum\_{i=2}^l \rho\_i \phi\_i(k)\Delta u(k-i+1)}{\lambda + \left|\hat{\phi}\_1(k)\right|^2} \\ &+ \frac{\phi\_1(k)\rho\_1(s^\*(k+1) - s(k) - \Delta\ell(k+1))}{\lambda + \left|\phi\_1(k)\right|^2}, \end{split} \tag{12}$$

*ϕ*ˆ1(*k*) = *ϕ*ˆ1(1), if *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*)- ≤ *ε*, or -Δ*Ul*(*k* − 1)-≤ *ε*, or sign(*ϕ*ˆ1(*k*)) = sign(*ϕ*ˆ1(1)), (13)

$$\mathfrak{e}(k) = \begin{cases} \ r\_{j\prime} \text{ if } \mathfrak{e}(k) > r\_{j\prime}s(k) = s\_{j} \\ -r\_{j\prime} \text{ if } \mathfrak{e}(k) \le -r\_{j\prime}s(k) = s\_{j}. \end{cases} \tag{14}$$

where *η* ∈ (0, 2], *μ* > 0, *γ* > 1, *λ* > 0, *ρ<sup>i</sup>* ∈ (0, 1], *i* ∈ {1, ··· , *l*}, *j* ∈ {1, ··· , *N*}; *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*) is the estimation vector of PG *ϕ<sup>f</sup>* ,*l*(*k*); *ε* denotes a small positive constant; *ϕ*ˆ1(1) is the initial value of *ϕ*ˆ1(*k*); the algorithm (13) is the reset algorithm of the PG estimation vector, and the algorithm (14) denotes the limitation mechanism of classification-metric deviation.

It is known from the above algorithms that the PG estimation vector directly affects the quality of the control scheme. In order to enhance the time-varying parameters' tracking ability for the PG estimation (10), it is necessary to add the reset algorithm (13). The limitation mechanism (14) is added to ensure that the deviation within one pattern class is not greater than the corresponding pattern class radius. The pseudo-order *l* is supposed to be less than or equal to the sum of the input and output orders (*m* + *n*). A large number of experiments show that the lower the system complexity, the smaller the value of *l* can be. On the contrary, the higher the system complexity, the greater the *l* should be. It is obvious that the proposed P-PFDL-MFAC algorithms in this work degenerate to the P-CFDL-MFAC algorithms designed in [30] when *l* = 1.

#### **4. Performance of the Closed-Loop System**

The focus of this section is to analyze the performance of the closed-loop tracking control system, i.e., to prove the tracking error bounded stability of the closed-loop control system. Before this, the following assumptions and lemmas are proposed.

**Assumption 4.** *Considering the nonlinear system (2), for any desired bounded output* ¯*dx*<sup>∗</sup> (*k* + 1)*, a bounded input u*∗(*k*) *always exists and it can make the system output equal to* ¯*dx*<sup>∗</sup> (*k* + 1)*.*

**Assumption 5.** *The signal of the first element of the PG vector ϕ<sup>f</sup>* ,*l*(*k*) *is assumed to be known and unchanged at any time k with* -Δ*Ul*(*k*)- = 0*, i.e., ϕ*1(*k*) ≥  > 0 *(or ϕ*1(*k*) ≤  < 0*), is a small positive constant. In this work, in order to simplify the derivation of the conclusion, it is always assumed that ϕ*1(*k*) ≥  > 0 *without loss of generality.*

**Lemma 2** ([22])**.** *Let*

*A* = ⎡ ⎢ ⎢ ⎢ ⎣ *a*<sup>1</sup> *a*<sup>2</sup> ··· *ap* 1 0 ··· 0 *... ... . . .* 1 0 ⎤ ⎥ ⎥ ⎥ ⎦ (*p*×*p*) .

*If* ∑*<sup>p</sup> <sup>i</sup>*=<sup>1</sup> |*ai*| < 1*, then s*(*A*) < 1*, where s*(*A*) *is the spectral radius of A.*

**Lemma 3** ([17])**.** *Let <sup>A</sup>* <sup>∈</sup> <sup>R</sup>*p*×*p. For any given <sup>ε</sup>* <sup>&</sup>gt; <sup>0</sup>*, there exists an induced consistent matrix norm such that* -*A<sup>v</sup>* ≤ *s*(*A*) + *ε*, *where s*(*A*) *has the same meaning as Lemma* 2*.*

It is known to all that Assumption 4 is a necessary condition for the design and solution of the control problem, and it also shows that the output of the system (2) is controllable. Many plants satisfy the condition of Assumption 5 to some extent, and its actual physical background is also very clear, i.e., the plant's output increasing or decreasing corresponds to the control input increasing or decreasing. Next, our main results will be proven.

**Lemma 4.** *For the system (2) and (3) using the P-PFDL-MFAC scheme (10)–(14) under Assumptions 2–5, ϕ*ˆ *<sup>f</sup>* ,*l*(*k*)*is bounded.*

**Proof of Lemma 4.** When -Δ*Ul*(*k* − 1)- ≤ *ε*, it is obvious that *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*) is bounded from the reset algorithm (13) of the P-PFDL-MFAC scheme. When -Δ*Ul*(*k* − 1)- > *ε*, subtracting *ϕ<sup>f</sup>* ,*l*(*k*) in both sides of Equation (10) obtains

$$\begin{split} \bar{\boldsymbol{\sigma}}\_{f,l}(k) &= \bar{\boldsymbol{\sigma}}\_{f,l}(k-1) - \boldsymbol{\sigma}\_{f,l}(k) + \boldsymbol{\sigma}\_{f,l}(k-1) + \frac{\eta \Delta l I\_l(k-1) \Delta s(k)}{\mu + \|\Delta l I\_l(k-1)\|^2} \\ &- \frac{\eta \Delta l I\_l(k-1) \boldsymbol{\Phi}\_{f,l}^T(k-1) \Delta l I\_l(k-1)}{\mu + \|\Delta l I\_l(k-1)\|^2} \\ &= \left[ I - \frac{\eta \Delta l I\_l(k-1) \Delta l I\_l^T(k-1)}{\mu + \|\Delta l I\_l(k-1)\|^2} \right] \bar{\boldsymbol{\sigma}}\_{f,l}(k-1) - \boldsymbol{\sigma}\_{f,l}(k) + \boldsymbol{\sigma}\_{f,l}(k-1) \\ &+ \frac{\eta \Delta l (k) \Delta l I\_l(k-1)}{\mu + \|\Delta l I\_l(k-1)\|^2} \end{split} \tag{15}$$

where *ϕ*˜ *<sup>f</sup>* ,*l*(*k*) = *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*) − *ϕ<sup>f</sup>* ,*l*(*k*).

Taking the norm on both sides of (15) and using Lemma 1, |*e*ˆ(*k*)| ≤ *rmax* yields

$$\|\|\vec{\varphi}\_{f,l}(k)\|\| \le 2b + 2\eta r\_{\text{max}} + \left\| \left[ I - \frac{\eta \Delta l I\_l(k-1) \Delta l I\_l^T(k-1)}{\mu + \|\Delta l I\_l(k-1)\|^2} \right] \vec{\varphi}\_{f,l}(k-1) \right\|.\tag{16}$$

Square the first term on the right of (16) and obtain the following inequality:

$$\begin{split} \left\lVert \left[ I - \frac{\eta \Delta l I\_l (k-1) \Delta l I\_l^T (k-1)}{\mu + \| \Delta l I\_l (k-1) \|^2} \right] \tilde{\boldsymbol{\sigma}}\_{f,l} (k-1) \right\rVert^2 &\leq \left\lVert \tilde{\boldsymbol{\sigma}}\_{f,l} (k-1) \right\rVert^2 + \\ \left( -2 + \frac{\eta \| \Delta l I\_l (k-1) \|^2}{\mu + \| \Delta l I\_l (k-1) \|^2} \right) \frac{\eta \left( \tilde{\boldsymbol{\sigma}}\_{f,l}^T (k-1) \Delta l I\_l (k-1) \right)^2}{\mu + \| \Delta l I\_l (k-1) \|^2} . \end{split} \tag{17}$$

Since *<sup>μ</sup>* <sup>&</sup>gt; 0 and *<sup>η</sup>* <sup>∈</sup> (0, 2], it can be obtained that <sup>−</sup><sup>2</sup> <sup>+</sup> *<sup>η</sup>*-Δ*Ul*(*k*−1)-2 *μ*+-Δ*Ul*(*k*−1)-<sup>2</sup> < 0, and it

is obvious that *<sup>η</sup> ϕ*˜*T f* ,*l* (*k*−1)Δ*Ul*(*k*−1) 2 *μ*+-Δ*Ul*(*k*−1)-<sup>2</sup> > 0. Thus, there must exist a constant 0 < *d*<sup>1</sup> < 1 that satisfies *<sup>I</sup>* <sup>−</sup> *<sup>η</sup>*Δ*Ul*(*k*−1)Δ*U<sup>T</sup> <sup>l</sup>* (*k*−1) *μ*+-Δ*Ul*(*k*−1)-2 *ϕ*˜ *<sup>f</sup>* ,*l*(*k* − 1) ≤ *d*<sup>1</sup> *ϕ*˜ *<sup>f</sup>* ,*l*(*<sup>k</sup>* <sup>−</sup> <sup>1</sup>) . It can be further deduced that

$$\begin{split} \|\bar{\varrho}\_{f,l}(k)\| &\leq d\_{1} \|\bar{\varrho}\_{f,l}(k-1)\| + 2b + 2\eta r\_{\max} \\ &\leq d\_{1}^{2} \|\bar{\varrho}\_{f,l}(k-1)\| + d\_{1}(2b + 2\eta r\_{\max}) + 2b + 2\eta r\_{\max} \\ &\leq \cdots \leq d\_{1}^{k-1} \|\bar{\varrho}\_{f,l}(1)\| + \frac{(2b + 2\eta r\_{\max})(1 - d\_{1}^{k-1})}{1 - d\_{1}}. \end{split} \tag{18}$$

In view of (18), *ϕ*˜ *<sup>f</sup>* ,*l*(*k*) is bounded, since *ϕ<sup>f</sup>* ,*l*(*k*) is bounded; thus, *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*) is bounded.

**Theorem 1.** *For system (2) and (3) using the P-PFDL-MFAC scheme (10)–(14) under Assumptions 3–6 with the desired signal s*∗(*k* + 1) = *s*∗ = *const, if the controller parameters meet the following conditions*

$$\text{(1)}\quad\text{letting }\overline{\rho}\_1 = \frac{\gamma^2 \rho\_1}{\gamma^2 - 1 + \rho\_1} \text{ and } \overline{\rho}\_1 \in (0, 1];$$

$$\text{(2)}\quad\text{letting }\mathfrak{p}\_{i} = \frac{(\gamma^{2}-1)\mathfrak{p}\_{i} + \mathfrak{p}\_{1}}{\gamma^{2}-1+\mathfrak{p}\_{1}}\text{ and }\mathfrak{p}\_{i} \in (0,1], i = 2, \cdots, l;$$

*(3) letting <sup>λ</sup>*¯ <sup>=</sup> (*γ*2−1)*<sup>λ</sup> <sup>γ</sup>*2−1+*ρ*<sup>1</sup> *, and there exists a λ*¯ *min such that λ*¯ > *λ*¯ *min,*

*then the closed-loop control system guarantees that*

$$\lim\_{k \to \infty} |s^\*-s(k+1)| \le M\_{\prime}$$

*where M is a constant and M* > 0*.*

**Proof of Theorem 1.** Substituting the classification-metric deviation estimation algorithm (11) into control algorithm (12), one has

$$u(k) = u(k-1) + \frac{\frac{\gamma^2 \rho\_1}{\gamma^2 - 1 + \rho\_1} \phi\_1(k)(\mathbf{s}^\* - \mathbf{s}(k))}{\frac{(\gamma^2 - 1)\lambda}{\gamma^2 - 1 + \rho\_1} + |\phi\_1(k)|^2} - \frac{\phi\_1(k) \sum\_{l=2}^l \frac{(\gamma^2 - 1)\rho\_l}{\gamma^2 - 1 + \rho\_1} \phi\_l(k) \Delta u(k - i + 1)}{\frac{(\gamma^2 - 1)\lambda}{\gamma^2 - 1 + \rho\_1} + |\phi\_1(k)|^2}. \tag{19}$$

Given *ρ*¯1, *ρ*¯*i*, *λ*¯ , Equation (19) can be written as

$$u(k) = u(k-1) + \frac{\bar{\rho}\_1 \hat{\phi}\_1(k)(s^\*-s(k))}{\bar{\lambda} + |\phi\_1(k)|^2} - \frac{\phi\_1(k) \sum\_{i=2}^l \bar{\rho}\_i \phi\_i(k) \Delta u(k-i+1)}{\bar{\lambda} + |\phi\_1(k)|^2},\tag{20}$$

where *ρ*¯*<sup>i</sup>* ∈ (0, 1], *i* = 1, ··· , *l*.

Since *<sup>γ</sup>* <sup>&</sup>gt; 1, *<sup>λ</sup>* <sup>&</sup>gt; 0 and *<sup>ρ</sup>*<sup>1</sup> <sup>∈</sup> (0, 1], thus *<sup>λ</sup>*¯ <sup>&</sup>gt; 0. It is known from Lemma 4 that *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*) is bounded and noted that *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*)- ≤ *b*1; here, *b*<sup>1</sup> is a positive constant. Given *ϕ*ˆ *<sup>f</sup>* ,*l*(*k*)- ≤ *b*<sup>1</sup> , *ϕ<sup>f</sup>* ,*l*(*k*)-<sup>≤</sup> *<sup>b</sup>*, *<sup>γ</sup>* <sup>&</sup>gt; 1, *<sup>λ</sup>* <sup>&</sup>gt; 0, *<sup>ρ</sup><sup>i</sup>* <sup>∈</sup> (0, 1], *<sup>ρ</sup>*¯*<sup>i</sup>* <sup>∈</sup> (0, 1], *<sup>λ</sup>*¯ <sup>&</sup>gt; 0, there exist bounded constants *Wi*, *i* ∈ {1, 2, 3, 4, 5} such that the following inequalities (21)–(25) hold when *λ*¯ > *λ*¯ *min*.

Letting *<sup>λ</sup>*¯ <sup>&</sup>gt; *<sup>λ</sup>*¯ *min* <sup>≥</sup> *<sup>b</sup>*<sup>2</sup> and using inequality *<sup>x</sup>*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*<sup>2</sup> <sup>≥</sup> <sup>2</sup>*xy* , one obtains

$$\left|\frac{\hat{\varphi}\_1(k)}{\overline{\lambda} + |\hat{\varphi}\_1(k)|^2}\right| \le \left|\frac{\hat{\varphi}\_1(k)}{2\sqrt{\lambda}|\hat{\varphi}\_1(k)|}\right| < \left|\frac{1}{2\sqrt{\lambda}\_{\text{min}}}\right| = W\_1 < \frac{0.5}{b},\tag{21}$$

$$0 < \mathcal{W}\_2 \le \left| \frac{\phi\_1(k)\phi\_i(k)}{\overline{\lambda} + |\phi\_1(k)|^2} \right| \le b \left| \frac{\phi\_1(k)}{2\sqrt{\mathcal{A}}|\phi\_1(k)|} \right| < 0.5,\tag{22}$$

$$\mathcal{W}\_1 \|\!\!\varphi\_{f,l}(k)\|\!\!\/=\mathcal{W}\_{\mathfrak{J}}<0.5.\tag{23}$$

From the inequalities (22) and (23), it is deduced that

$$\mathcal{W}\_2 + \mathcal{W}\_3 < \mathbf{1}.\tag{24}$$

Letting {∑*<sup>l</sup> i*=2 *ϕ*ˆ1(*k*)*ϕ*ˆ*i*(*k*) *<sup>λ</sup>*¯ <sup>+</sup>|*ϕ*ˆ1(*k*)|<sup>2</sup> } 1 *<sup>l</sup>*−<sup>1</sup> ≤ *<sup>W</sup>*<sup>4</sup> and choosing *<sup>ρ</sup>*¯*max* = max*i*=1,··· ,*<sup>l</sup> <sup>ρ</sup>*¯*i*, one has

$$\sum\_{i=2}^{I} \overline{\rho}\_{i} \left| \frac{\hat{\varrho}\_{1}(k)\hat{\varrho}\_{i}(k)}{\overline{\lambda} + |\phi\_{1}(k)|^{2}} \right| \leq \overline{\rho}\_{\text{max}} \sum\_{i=2}^{I} \left| \frac{\hat{\varrho}\_{1}(k)\hat{\varrho}\_{i}(k)}{\overline{\lambda} + |\phi\_{1}(k)|^{2}} \right| \leq \overline{\rho}\_{\text{max}} W\_{4}^{l-1} = W\_{5} < 1. \tag{25}$$

Defining tracking error *w*(*k*) = *s*<sup>∗</sup> − *s*(*k*) and letting

$$A(k) = \begin{bmatrix} -\frac{\rho\_2 \Phi\_1(k)\Phi\_2(k)}{\bar{\lambda} + |\bar{\varphi}\_1(k)|^2} & -\frac{\rho\_3 \Phi\_1(k)\Phi\_3(k)}{\bar{\lambda} + |\bar{\varphi}\_1(k)|^2} & \cdots & -\frac{\rho\_l \Phi\_1(k)\Phi\_l(k)}{\bar{\lambda} + |\bar{\varphi}\_1(k)|^2} & 0\\ 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix},\tag{26}$$

the control algorithm (12) can be written as

$$\begin{split} \Delta \mathcal{U}\_l(k) &= [\Delta u(k), \cdot, \cdot, \Delta u(k - l + 1)]^T \\ &= A(k) \left[ \Delta u(k - 1), \cdot, \cdot, \Delta u(k - l) \right]^T + \frac{\bar{\rho}\_1 \phi\_1(k)}{\bar{\lambda} + |\phi\_1(k)|^2} \mathcal{C}w(k), \end{split} \tag{27}$$

where *C* = [1, 0, ··· , 0] *<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*<sup>l</sup>* . The secular equation of *A*(*k*) is

$$z^I + \frac{\bar{\rho}\_2 \phi\_1(k) \phi\_2(k)}{\bar{\lambda} + |\hat{\phi}\_1(k)|^2} z^{I-1} + \dots + \frac{\bar{\rho}\_l \phi\_1(k) \phi\_l(k)}{\bar{\lambda} + |\hat{\phi}\_1(k)|^2} z = 0.$$

From Lemma 2 and inequality (25), one has |*z*| < 1 and obtains

$$|z|^{l-1} \le \sum\_{i=2}^{l} \bar{\rho}\_i \left| \frac{\phi\_1(k)\phi\_i(k)}{\bar{\lambda} + |\phi\_1(k)|^2} \right| \le \bar{\rho}\_{\max} \mathcal{W}\_4^{l-1} < 1. \bar{\rho}\_i$$

Further, it can be deduced that |*z*| ≤ *ρ*¯ 1 *<sup>l</sup>*−<sup>1</sup> *maxW*4. From Lemma 3, one can obtain -*A*(*k*)*<sup>v</sup>* ≤ *s*(*A*(*k*)) + *ε* ≤ *ρ*¯ 1 *<sup>l</sup>*−<sup>1</sup> *maxW*<sup>4</sup> < 1. According to the definition of *Ul*(*k*), it is clear that Δ*Ul*(0) = 0. Letting *d*<sup>2</sup> = *ρ*¯ 1 *<sup>l</sup>*−<sup>1</sup> *maxW*<sup>4</sup> and taking the norm on both sides of (27), one obtains

$$\begin{split} \|\Delta \mathcal{U}\_{l}(k)\| &\leq \|A(k)\|\_{\mathbb{P}} \|\Delta \mathcal{U}\_{l}(k-1)\| + \bar{\rho}\_{1} \left| \frac{\phi\_{1}(k)}{\bar{\lambda} + |\phi\_{1}(k)|^{2}} \right| |w(k)| \\ &\leq d\_{2} \|\Delta \mathcal{U}\_{l}(k-1)\| + \bar{\rho}\_{1} \mathcal{W}\_{l} |w(k)| \leq \cdot \cdot = \bar{\rho}\_{1} \mathcal{W}\_{1} \sum\_{i=1}^{k} d\_{2}^{k-i} |w(k)|. \end{split} \tag{28}$$

From Lemma 1 and Equation (27), one has

$$\begin{split} w(k+1) &= s^\* - s(k+1) = s^\* - \bar{d}\mathbf{x}(k+1) - e(k+1) \\ &= w(k) - \Delta e(k+1) - q^T\_{f,l}(k)\Delta l I\_l(k) \\ &= \left[ 1 - \frac{\bar{\rho}\_1 \hat{\rho}\_1(k)\varrho\_1(k)}{\bar{\lambda} + |\Phi\_1(k)|^2} \right] w(k) - q^T\_{f,l}(k)A(k)\Delta l I\_l(k-1) - \Delta e(k+1). \end{split} \tag{29}$$

Choosing a reasonable *ρ*¯1, one can obtain

$$\left|1-\frac{\bar{\rho}\_1\phi\_1(k)\,\rho\_1(k)}{\bar{\lambda}+|\hat{\phi}\_1(k)|^2}\right| = \left|1-\left|\frac{\bar{\rho}\_1\phi\_1(k)\,\rho\_1(k)}{\bar{\lambda}+|\hat{\phi}\_1(k)|^2}\right|| \le 1-\bar{\rho}\_1\mathcal{W}\_2 = d\eta < 1. \tag{30}$$

From the above inequality and |*e*(*k*)| ≤ *rmax*, taking the norm on both sides of the Equation (29), one obtains

$$\begin{split} |w(k+1)| &< d\_3|w(k)| + d\_2 ||\varphi\_{f,l}(k)|| \|\Delta l\_l(k-1)|| + 2r\_{\max} < \cdots \\ &< d\_3^k |w(1)| + d\_2 \sum\_{i=1}^{k-1} d\_3^{k-1-i} ||\varphi\_{f,l}(i+1)|| \|\Delta l\_l(i)\| + 2r\_{\max} \sum\_{i=1}^{k-1} d\_3^{k-1-i} \\ &< d\_3^k |w(1)| + 2r\_{\max} \sum\_{i=1}^{k-1} d\_3^{k-1-i} + d\_2 \sum\_{i=1}^{k-1} d\_3^{k-1-i} ||\varphi\_{f,l}(i+1)|| |\rho\_1 \mathcal{W}\_1 \sum\_{j=1}^i d\_2^{i-j} |w(j)|. \end{split} \tag{31}$$

Letting *d*<sup>4</sup> = *ρ*¯1*W*3, it is clear that *d*<sup>4</sup> < 1. The inequality (31) can be recorded as

$$|w(k+1)| < d\_3^k |w(1)| + d\_2 d\_4 \sum\_{i=1}^{k-1} d\_3^{k-1-i} \sum\_{j=1}^i d\_2^{i-j} |w(j)| + \frac{2r\_{\max}(1 - d\_3^{k-1})}{1 - d\_3}.\tag{32}$$

Letting

$$\mathcal{g}(k+1) = d\_3^k |w(1)| + d\_2 d\_4 \sum\_{i=1}^{k-1} d\_3^{k-1-i} \sum\_{j=1}^i d\_2^{i-j} |w(j)|.$$

it is obvious that *g*(2) = *d*3|*w*(1)|. One can see that if *g*(*k* + 1) is bounded, then *w*(*k*) is bounded.

Next, the boundedness of *g*(*k* + 1) will be proven.

$$\begin{split} g(k+2) &= d\_3^{k+1}|w(1)| + d\_2 d\_4 \sum\_{i=1}^k d\_3^{k-i} \sum\_{j=1}^i d\_2^{j-j} |w(j)| \\ &= d\_3 \mathfrak{g}(k+1) + d\_4 d\_2^k |w(1)| + \dots + d\_4 d\_2^2 |w(k-1)| + d\_4 d\_2 |w(k)| \\ &< d\_3 \mathfrak{g}(k+1) + d\_4 d\_2^k |w(1)| + \dots + d\_4 d\_2 \mathfrak{g}(k) \\ &+ d\_4 d\_2^2 |w(k-1)| + d\_4 d\_2 \frac{2r\_{\max}(1 - d\_3^{k-2})}{1 - d\_3}. \end{split} \tag{33}$$

$$\begin{array}{c} \text{Note that } \bar{h}(k) = d\_3 \mathfrak{g}(k+1) + d\_4 d\_2^k |w(1)| + \dots + d\_4 d\_2^2 |w(k-1)| + d\_4 d\_2 \mathfrak{g}(k). \text{ Since } \bar{h} = 1 - \bar{\rho}\_1 \mathbb{W}\_2 > \bar{\rho}\_1 (\mathbb{W}\_2 + \mathbb{W}\_3) - \bar{\rho}\_1 \mathbb{W}\_2 = \bar{\rho}\_1 \mathbb{W}\_3 = d\_4 \text{, one obtains} \end{array}$$

$$\begin{split} &h(k) < d\_3 \mathfrak{z}(k+1) + d\_4 d\_2^k |w(1)| + \dots + d\_4 d\_2^2 |w(k-1)| + d\_3 d\_2 \mathfrak{z}(k) \\ &< d\_3 \mathfrak{z}(k+1) + d\_4 d\_2^k |w(1)| + \dots + d\_4 d\_2^2 |w(k-1)| \\ &+ d\_3 d\_2 \left[ d\_3^{k-1} |w(1)| + d\_2 d\_4 \sum\_{i=1}^{k-2} d\_3^{k-2-i} \sum\_{j=1}^i d\_2^{j-j} |w(j)| \right] \\ &= d\_2 \mathfrak{z}(k+1). \end{split} \tag{34}$$

From the inequalities (33) and (34), one has

$$
\lg(k+2) \le (d\_2 + d\_3)\lg(k+1) + d\_4 d\_2 \frac{2r\_{\max}(1 - d\_3^{k-2})}{1 - d\_3}.
$$

Since *d*<sup>2</sup> + *d*<sup>3</sup> = 1 − *ρ*¯1*W*<sup>2</sup> + *ρ*¯ 1 *<sup>l</sup>*−<sup>1</sup> *maxW*4, by choosing the reasonable *ρ*¯*i*, *i* = 1, ··· , *l*, it exits *d*<sup>2</sup> + *d*<sup>3</sup> = *d*<sup>5</sup> ∈ (0, 1) and one obtains

$$d\_1 g(k+2) \le d\_5 g(k+1) + d\_4 d\_2 \frac{2r\_{\text{max}}}{1 - d\_3} \le \dots \le d\_5^k g(2) + d\_4 d\_2 \frac{2r\_{\text{max}}}{1 - d\_3} \frac{1 - d\_5^k}{1 - d\_5}.\tag{35}$$

.

It is clear that *g*(*k*) is bounded convergent; thus, the tracking error *w*(*k*) is bounded convergent, i.e., lim*k*→<sup>∞</sup> |*w*(*k*)| ≤ *<sup>M</sup>*, *<sup>M</sup>* is a positive constant.

**Remark 3.** *The contraction mapping principle is utilized to prove the bounded convergence in this work, and many inequalities are employed to handle the mapping relationships in Lemma 4 and Theorem 1. A critical technique is to let λ, γ, and ρ<sup>i</sup> take reasonable values that can guarantee the existence of constants W*1*, W*2*, W*3*, W*4*, W*5*, λ*¯ *, γ*¯*, ρ*¯*i, d*1*, d*2*, d*3*, d*4*, and d*<sup>5</sup> *to make the inequalities used in the above derivations hold.*

**Remark 4.** *It is obvious that the desired tracking target is an arbitrary bounded constant s*∗ *in Theorem 1. In fact, for the closed-loop control system based on pattern moving, the desired tracking target should be one or some specific pattern classes (dxi), i.e., one or some specific pattern class centers (s*<sup>∗</sup> = *si, i* = 1, ··· , *N). Therefore, instead of focusing on each specific value of the system output, the P-PFDL-MFAC method focuses on whether the system outputs belong to one or some specific pattern classes, and this is the most significant difference between the method designed in this work and the model free adaptive quantization control method proposed in [35,36]. From this point of view, under the control input and output disturbance, even if the implicit metric value of the pattern class to which the system outputs belong satisfies* <sup>|</sup> ¯*dx*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> *<sup>s</sup>*∗| ≤ *ri when the desired target s*<sup>∗</sup> = *si, it is still considered that the system's tracking error is zero.*

**Remark 5.** *The designed P-PFDL-MFAC method is employed for the considered system (2) and (3), which corresponds to a practical SIMO system (1). When the system is under the control input u*(*k*) *at time instant k, the output vector* [*y*1(*k* + 1), ··· , *yl*(*k* + 1)] *is obtained, and then s*(*k* + 1) *is obtained by feature extraction T*(·)*, pattern classification M*(·)*, the class center explicit metric D*(·) *with the real-time output data* [*y*1(*k* + 1), ··· , *yl*(*k* + 1)]*, and a large amount of offline historical data. Generally speaking, the P-PFDL-MFAC method can be considered a novel data-driven method based on offline historical data and online real-time data, and this is a major difference from the traditional MFAC methods.*

#### **5. Simulation**

Two examples are given to demonstrate the feasibility and effectiveness of the achieved algorithms in this section. In the simulation example of reference [37], the speed control of a Stanford manipulator's joint 4 proposed in [38] was discussed. It considered that the controlled object is a discrete-time system with jump parameters while the load changes. In the first example below, this discrete-time system is also taken as the consideration object, and the designed P-PFDL-MFAC scheme is implemented. Example 2 is a SIMO nonlinear discrete-time numerical case. In this simulation case, the designed control scheme is adopted, and the control effects with different pseudo-orders are compared.

**Example 1.** *Consider a SISO discrete-time system with jump parameters*

$$y(k) = a\_2(k)y(k-2) + b\_0(k)
\mu(k-1) + b\_1(k)
\mu(k-2) + g(k) + e(k),\tag{36}$$

*where y*(*k*) *is the system output, which denotes the speed of a Stanford manipulator's joint* 4*; u*(*k*) *is the system input, which denotes the motor's voltage and satisfies u*(*k*) ∈ [0, 10]*; e*(*t*) *denotes the system random noise and it satisfies that* |*e*(*k*)| ≤ 0.01*; g*(*k*) *is considered as a constant and g*(*t*) = 0.25*; b*1(*k*) *is also a constant and b*1(*k*) = 0.2*; the other two system jump parameters are as follows:*

$$a\_2(k) = \begin{cases} -0.9, k \le 200; \\ -0.75, 200 < k \le 400; \\ -0.9, 400 < k \le 600, \end{cases}$$

*and*

$$b\_0(k) = \begin{cases} 0.4, k \le 200; \\ 0.35, 200 < k \le 400; \\ 0.4, 400 < k \le 600. \end{cases}$$

*The control goal of our designed scheme is that the outputs belong to one or some special pattern classes, which is the most significant difference from the simulation in [37]. Firstly, a large number of outputs obtained under effective control inputs are divided into several pattern classes. Then, one or some desired pattern classes are taken as the targets of system control.*

Step 1: Classification (*M*(·)) and metrics (*D*(·), *<sup>D</sup>*¯ (·)) of massive offline data. Here, 600 evenly distributed inputs are taken and the corresponding outputs are obtained. A modified quantized control classification and class center explicit metric method (*M*(·), *D*(·)) [34] is adopted and described as follows.

$$s(k) = D(M(y(k))) = \begin{cases} y\_0(k), & \text{if } T1\_i < y(k) \le T2\_{i\prime} \\ 0, & \text{if } -T\_N < y(k) \le T\_{N\prime} \\ -y\_0(k), \text{if } -T2\_i < y(k) \le -T1\_{i\prime} \end{cases} \tag{37}$$

where *T*1*<sup>i</sup>* = <sup>1</sup> <sup>1</sup>+*Δκi*; *<sup>T</sup>*2*<sup>i</sup>* <sup>=</sup> <sup>1</sup> <sup>1</sup>−*Δκi*; *TN* <sup>=</sup> <sup>1</sup> 1+*Δρ<sup>N</sup>* <sup>0</sup> *<sup>κ</sup>*0; *<sup>y</sup>*0(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) = <sup>1</sup>+*ρ*<sup>0</sup> <sup>4</sup> *<sup>κ</sup>i*(*ρi*−<sup>1</sup> <sup>0</sup> <sup>+</sup> *<sup>ρ</sup><sup>i</sup>* <sup>0</sup>); *<sup>Δ</sup>* <sup>=</sup> <sup>1</sup>−*ρ*<sup>0</sup> 1+*ρ*<sup>0</sup> ; *κ<sup>i</sup>* = *ρ<sup>i</sup>* <sup>0</sup>*κ*0; *ρ*<sup>0</sup> ∈ (0, 1); *κ*<sup>0</sup> is the maximum working range of *y*(*k*) (*κ*<sup>0</sup> ≥ max{|*y*(*k*)|}); *N* denotes the number of pattern classes; *i* = 1, 2, ··· , *N* − 1.

Given the upper limit of the initial class radius *r*<sup>0</sup> at the working point 0 and other parameters such as *ρ*<sup>0</sup> and *κ*0, one can obtain *L* ≥ ln(*r*<sup>0</sup> (1+*Δ*) *<sup>κ</sup>*<sup>0</sup> ) ln *<sup>ρ</sup>*<sup>0</sup> , and the output sequence {*y*(*k*)} is divided into 2*<sup>L</sup>* <sup>+</sup> 1 segments. Furthermore, *<sup>N</sup>* <sup>=</sup> <sup>2</sup>*<sup>L</sup>* <sup>+</sup> 1, *si*, *ri* <sup>=</sup> <sup>1</sup>+*ρ*<sup>2</sup> <sup>4</sup>*<sup>ρ</sup>* and class threshold *Ci* can be obtained, respectively, *i* = 1, ··· , *N*. The parameter settings of the adopted classification method are *ρ*<sup>0</sup> = 0.4, *κ*<sup>0</sup> = 15, *r*<sup>0</sup> = 0.2. The distribution curves of {*u*(*k*)}, {*y*(*k*)}, and {*s*(*k*)} are shown in Figure 1. Table 1 shows the property values of each pattern class.


**Table 1.** Property values of pattern class.

**Figure 1.** The curves of I/O data and class centers.

**Remark 6.** *To the best of our knowledge, there are many clustering and classification algorithms in statistical pattern recognition, such as ISODATA, K-means, C-means, and so on. A class center explicit metric and modified quantized control classification method is adopted in this work. As mentioned in [2], the product quality is directly related to the working conditions. Therefore, the parameter settings of condition classification are determined by the result of product quality clustering. Here, it is assumed that the first principal component information y*(*k*) ∈ (0.2688, 0.6720] *corresponds to good product quality, so the initial parameters (ρ*<sup>0</sup> = 0.4*, κ*<sup>0</sup> = 15*, r*<sup>0</sup> = 0.2*) are configured to ensure that the working condition data y*(*k*) ∈ (0.2688, 0.6720] *belong to one pattern class.*

Step 2: A pattern-moving-based system dynamics description is established with the obtained property values and data sets {*u*(*k*)}, ¯*dx*(*k*) , and {*s*(*k*)}.

$$\begin{cases} \bar{d}\mathbf{x}(k) = f(\bar{d}\mathbf{x}(k-1), \dots, \bar{d}\mathbf{x}(k-n\_y), u(k-1), \dots, u(k-n\_u)), \\\\ s(k) = D(M(\bar{d}\mathbf{x}(k)) = \begin{cases} -7.3500, \bar{d}\mathbf{x}(k) \in (-10.5, -4.2], \\ \vdots \\ 0.0000, \bar{d}\mathbf{x}(k) \in (-0.1075, 0.1075], \\ \vdots \\ \end{cases} \end{cases} \tag{38}$$

where *f*(·) is an unknown nonlinear system function; *nu*, *ny* denote the unkown input and output orders of *f*(·), respectively.

Step 3: Application of the control scheme. Nine pattern classes are obtained and the designed P-PFDL-PMFAC scheme (10)–(14) is employed to track the following targets.

$$s^\*(k) = 0.4704\_\*$$

where *s*∗ = 0.4704 denotes that the object is pattern class 8.

Set the initial conditions as *y*(1 : 2) = 0, *e*(1 : 2) = 0, *u*(1 : 2) = 0, *ϕ*ˆ1(2) = 1, *ϕ*ˆ2(1 : 2) = 0, *ε* = 10<sup>−</sup>5, *s*(1 : 2) = 0. The controller parameters are set as *γ* = 10, *λ* = 0.01, *μ* = 1, *η* = 0.5, *ρ*<sup>1</sup> = *ρ*<sup>2</sup> = 0.5, *l* = 2 and the resetting initial value is *ϕ*ˆ1(1) = 0.5. Figure 2 shows the system output process, and Figure 3 shows the curves of control input, PG estimation values, and deviation. From the controlled output of the system, it can be seen that although it has undergone drastic adjustment at the beginning, it can track the target quickly and achieve a good tracking effect.

**Figure 2.** The curves of desired class center, original output, and its corresponding class center.

**Figure 3.** The curves of control input, PG estimation values, and deviation.

**Example 2.** *A single input and three outputs of the nonlinear discrete-time system are given as follows.*

$$\begin{cases} y\_1(k+1) &= 1.2\sin(0.5y\_1(k)) + u^2(k) + \frac{u(k)}{1+u^2(k)} + u(k-1) + d(k), \\ y\_2(k+1) &= 1.3\sin(0.5y\_2(k)) + 0.2y\_2(k-1) + \frac{u(k)}{1+u^2(k)} + 0.5u(k-1) + d(k), \\ y\_3(k+1) &= 1.4\sin(0.5y\_3(k)) + 0.5u^2(k) + \frac{u(k)}{1+u^2(k)} + u(k-1) + d(k), \end{cases} \tag{39}$$

*where yi*(*k*) *denotes one of the three outputs, i* = 1, 2, 3*; d*(*k*) *is the Gaussian white noise and <sup>d</sup>*(*k*) ∼ N (0, 0.012)*; <sup>u</sup>*(*k*) *denotes the system input and <sup>u</sup>*(*k*) ∈ [−2, 2]*; the system is merely employed to produce the I/O data with unknown system structure, orders, and parameters.*

Feature extraction (*T*(·)), classification (*M*(·)), and metrics (*D*(·), *<sup>D</sup>*¯ (·)) of massive offline data. Here, 1000 evenly distributed inputs are taken and the corresponding outputs are obtained. The outputs are normalized and the PCA technology is employed to deal with them. One can obtain the first principal component information {*y*(*k*)} (the contribution rate: 85.4518% > 85%). The same classification-metrics method (37) as in Example 1 is adopted. The parameter settings of the adopted classification method are *ρ*<sup>0</sup> = 0.4, *κ*<sup>0</sup> = 5, *r*<sup>0</sup> = 0.2. The distribution curves of {*u*(*k*)}, {*yi*(*k*)}, {*y*(*k*)}, and {*s*(*k*)} are shown in Figure 4, *i* = 1, 2, 3. Table 2 shows the property values of each pattern class.

**Figure 4.** The curves of I/O data, PCA information, and class center.


**Table 2.** Property values of pattern class.

A pattern-moving-based system dynamics description is established as follows.

$$\begin{cases} \bar{d}\mathbf{x}(k+1) = f(\bar{d}\mathbf{x}(k), \dots, \bar{d}\mathbf{x}(k-n\_y), u(k), \dots, u(k-n\_u)), \\\\ \mathbf{s}(k+1) = D(M(\bar{d}\mathbf{x}(k+1))) = \begin{cases} -2.4500, \bar{d}\mathbf{x}(k+1) \in (-3.5, -1.4], \\ \vdots \\ 0.0000, \bar{d}\mathbf{x}(k+1) \in (-0.0896, 0.0896], \\ \vdots \\ \vdots \\ 2.4500, \bar{d}\mathbf{x}(k+1) \in (1.4, 3.5], \end{cases} \end{cases} \tag{40}$$

Nine pattern classes are obtained and the designed P-PFDL-PMFAC scheme (10)–(14) is employed to track the following targets.

$$\mathbf{s}^\*(k) = \begin{cases} 0.000, & 0 < k \le 500; \\ 0.980, 500 < k \le 1000. \end{cases}$$

where *s*∗ = 0, *s*∗ = 0.980 denote that the object is pattern class 5 and 8, respectively.

Set the initial conditions as *y*1(1 : 4) = 0, *y*2(1 : 4) = 0, *y*3(1 : 4) = 0, *e*(1 : 4) = 0, *u*(1 : 4) = 0, *ϕ*ˆ1(2 : 4) = 1, *ϕ*ˆ2(1 : 4) = 0, *ϕ*ˆ3(1 : 4) = 0, *ε* = 10−5, *s*(1 : 4) = 0 . The controller parameters are set as *γ* = 10, *λ* = 0.01, *μ* = 1, *η* = 0.5, *ρ*<sup>1</sup> = *ρ*<sup>2</sup> = *ρ*<sup>3</sup> = 0.5 and the resetting initial value is *ϕ*ˆ1(1) = 0.5. Figures 5–7 correspond to the curves of system input, outputs, PG estimation values, and deviation when the pseudo-order *l* is 1, 2, and 3, respectively. When *l* = 1, the P-PFDL-PMFAC scheme degenerates to the P-CFDL-MFAC method designed in [30], and the PG vector becomes a PPD. All three figures show that the target trajectory *s*∗(*k*) = 0.980 is well tracked. However, Figure 5 shows that the tracking effect of target trajectory *s*∗(*k*) = 0 is poor. Figure 6 shows that the tracking effect of target trajectory *s*∗(*k*) = 0 is slightly better, but there are also many cases where the tracking can not be achieved. It can be seen from Figure 7 that the target object *s*∗(*k*) = 0 is well tracked. The simulation results confirm that the value of pseudo-order should correspond to the complexity of the system, and they show that a reasonable pseudo-order can improve the control effect of the system. This numerical example illustrates that the designed scheme is a very feasible method for a class of nonlinear discrete-time systems when the outputs only need to be controlled to one or some specific pattern classes.

**Figure 5.** The curves of PPD estimation value, control input, deviation, and outputs with *l* = 1.

**Figure 6.** The curves of PG estimation values, control input, deviation, and outputs with *l* = 2.

**Figure 7.** The curves of PG estimation values, control input, deviation, and outputs with *l* = 3.

#### **6. Conclusions**

A novel P-PFDL-MFAC scheme is proposed by combining the pattern-moving-based system dynamics description with the traditional PFDL-MFAC approach for a class of unknown practical SIMO nonaffine nonlinear discrete-time systems. Obviously, this scheme can also be applied to nonlinear or linear time-varying SISO systems, as long as the purpose of system control is to make all outputs belong to one or some pattern classes. Due to the existence of classification-metric deviation, an improved cost function for a deviation estimation algorithm and an adaptive tracking control law is designed based on the saddle point theory of TP-ZSG. The bounded convergence of the closed-loop system's tracking error has been proven and the effectiveness of the P-PFDL-MFAC scheme has been validated via two simulation examples.

Although it can be seen from the simulation results that the control strategy proposed in this work has a good effect on the output disturbance, the robustness of data-driven control should also include the ability to deal with data dropout, which may be caused by sensor fault, transmission network failure, or actuator damage. Therefore, the next topic that needs to be focused on is the robustness of pattern-moving-based model free adaptive control in the case of missing data.

**Author Contributions:** Conceptualization, X.L. and Z.X.; methodology, Z.X.; software, X.L.; validation, X.L.; formal analysis, X.L.; investigation, X.L.; resources, X.L.; data curation, Z.X.; writing original draft preparation, X.L.; writing—review and editing, X.L.; visualization, X.L.; supervision, Z.X.; project administration, X.L.; funding acquisition, X.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (62076025).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Bin Kou 1,2, Shijie Guo 1,2,\* and Dongcheng Ren 1,2**


**Abstract:** Identifying the kinetic parameters of an industrial robot is the basis for designing a controller for it. To solve the problems of the poor accuracy and easy premature convergence of common bionic algorithms for identifying the dynamic parameters of such robots, this study proposed simulated annealing with similar exponential changes based on the beetle swarm optimization (SEDSABSO) algorithm. Expressions for the dynamics of the industrial robot were first obtained through the SymPyBotics toolkit in Python, and the required trajectories of excitation were then designed to identify its dynamic parameters. Following this, the search pattern of the global optimal solution for the beetle swarm optimization algorithm was improved in the context of solving for these parameters. The global convergence of the algorithm was improved by improving the iterative form of the number N of skinks in it by considering random perturbations and the simulated annealing algorithm, whereas its accuracy of convergence was improved through the class exponential change model. The improved beetle swarm optimization algorithm was used to identify the kinetic parameters of the Zhichang Kawasaki RS010N industrial robot. The results of experiments showed that the proposed algorithm was fast and highly accurate in identifying the kinetic parameters of the industrial robot.

**Keywords:** industrial robot; kinetic parameter identification; beetle swarm optimization algorithm; stochastic perturbation

#### **1. Introduction**

Kinetic parameters are the main factor influencing the control of fast and highly precise movements of industrial robots [1]. The process of identifying of their kinetic parameters is usually divided into a number of steps, such as kinetic modeling, designing the excitation trajectory, data acquisition, and identifying and verifying the kinetic parameters [2]. Gautier et al. used a two degree-of-freedom (DOF) robot as the object of study and applied the extended Kalman filter and the least-squares method to identify its parameters. Memar used the SCHUNK Powerball LWA 4P as an experimental object, constructed a dynamics model for it, and implemented the least-squares method to identify the dynamic parameters of the industrial robot [3]. However, the least-squares method is susceptible to measurement noise that lowers its accuracy [4]. Fu et al. [5] used the particle swarm optimization (PSO) algorithm with least squares to identify the kinetic parameters of a seven-DOF collaborative robot in Xinsong, but PSO can easily fall into the local optimum owing to poor population diversity in the late stage of processing that reduces the accuracy of identification of the parameters. Ding et al. [6] identified the dynamics of the robot by using the genetic algorithm, but the process of coding of the algorithm is cumbersome.

Summarizing the existing research, it is found that the identification of robot dynamics parameters is a high-dimensional function problem. Common algorithms are either the PSO accuracy is not high enough, or the genetic algorithm coding is more complicated, so a simple programming, strong anti-interference ability, and convergence accuracy are needed. High algorithm. BSO is a bionic algorithm recently proposed, which mainly uses the principles of PSO and the Beetle Antenna Search (BAS) algorithm, so it has the

**Citation:** Kou, B.; Guo, S.; Ren, D. A New Method for Identifying Kinetic Parameters of Industrial Robots. *Actuators* **2022**, *11*, 2. https:// doi.org/10.3390/act11010002

Academic Editors: Marco Carricato and Edoardo Idà

Received: 29 November 2021 Accepted: 21 December 2021 Published: 23 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

advantage of simple programming [7–9]. This article proposes SEDSABSO on the basis of BSO.SEDSABSO combined random perturbation-based behavior and the simulated annealing (SA) algorithm with the global optimal solution of the BSO [10,11] to improve its ability to search for global optimal particles as well as the manner of changes in N particles through an exponential decay model. This improved the convergence of BSO without affecting its computational complexity. Therefore, using the SEDSABSO algorithm for robot dynamic parameter identification will have the advantages of simple programming, high convergence accuracy, and fast iteration speed.

The rest of the paper is organized as follows. Section 2 describes the D-H and the dynamic parameters of the RS010N robot and then linearized them through the SymPyBotics toolkit in Python. In Section 3, the genetic algorithm toolbox in MATLAB was used to design the excitation trajectories required to identify the dynamic parameters of the robot. Section 4 describes the principle of the SEDSABSO algorithm and its process for the identification of robot dynamics parameters. In Section 5, the minimum set of parameters for the dynamics of the RS010N robot was identified and was used to compare the per-formance of the SEDSABSO, BSO, and LDWPSO algorithms. The results verified the accuracy and effectiveness of the proposed algorithm [12].

#### **2. Robotic arm Dynamics Model**

#### *2.1. RS010N Industrial Robots*

In this paper, we used the Zhichang Kawasaki RS010N robot as the research object. Figure 1 shows its structural configuration, and Figure 2 shows the configuration of its DH coordinates. The RS010N robot is a typical six-DOF industrial robot. Because the kinetic parameters of its three rear joints are much smaller than those of its three front joints, the former have a smaller influence on the accuracy of control of the robot's motion. Thus, we considered only the first three joints of the RS010N industrial robot here [13].

**Figure 1.** RS010N robot.

**Figure 2.** RS010N Robot DH coordinate system.

#### *2.2. Dynamics Modeling*

For an *n*-link robotic arm, the dynamics as modeled by the Newton-Euler method can be described as

$$\tau = D(q) \cdot \ddot{q} + \mathcal{C}(q, \dot{q}) + G(q) \tag{1}$$

In Equation (1), *τ* denotes the driving moment, *q* denotes the vector of the position of the joints, and . *<sup>q</sup>* and .. *q* denote vectors of the velocity and acceleration of the joints, respectively. *D*(*q*) denotes the inertial matrix, *C*(*q*, . *q*) denotes the Koch and centrifugal force terms, respectively, and *G*(*q*) denotes the gravitational force term [14].

According to the improved Newtonian–Eulerian dynamics, the above equation can be transformed into

$$
\pi = \Phi(q, \dot{q}, \ddot{q}) P \tag{2}
$$

where *Φ* is the observation matrix with a size of *n* ∗ 12*n*, and *P* denotes the vector of the inertial parameters of the robot arm.

#### *2.3. Minimum Set of Parameters for the Dynamical Model*

The parameters of the first three joint linkages of the RS010N six-freedom robot are shown in Table 1 below. *a*<sup>0</sup> = 350 mm, *a*<sup>1</sup> = 1035 mm.

**Table 1.** RS010N connecting rod parameters.


The kinetic parameters are shown in Table 2.


**Table 2.** Actual kinetic parameters.

We obtained the minimum set of parameters and related expressions for them through the SymPyBotics toolkit in Python.


(3)

According to Equation (3) above, the 36 kinetic parameters of the robot were linearized to 21, and the kinetic model can be written in the following form

$$
\pi = \Phi\_b(q, \dot{q}, \ddot{q}) p\_b \tag{4}
$$

where *Φ<sup>b</sup>* is the full-rank observation matrix and *Pb* denotes the underlying parameter vector of the robotic arm.

#### **3. Incentive Track Design**

The commonly used excitation trajectory when identifying the kinetic parameters of industrial robots is the finite-term Fourier series [15,16]. It can be expressed in the following form:

$$\ddot{\theta}\_i(t) = -\sum\_{n=1}^{N} a\_n n w\_f \sin(n w\_f t) - b\_n n w\_f \cos(n w\_f t) \tag{5}$$

$$\dot{\theta}\_i(t) = \sum\_{n=1}^{N} a\_n \cos(nw\_f t) - b\_n \sin(nw\_f t) \tag{6}$$

$$\theta\_i(t) = \theta\_0 + \sum\_{n=1}^{N} \frac{a\_n}{n w\_f} \sin(n w\_f t) - \frac{b\_n}{n w\_f} \cos(n w\_f t) \tag{7}$$

*wf* denotes the fundamental frequency of the Fourier series. Each Fourier series contains *an*, *bn*, and *θ*0. *N* represents the number of harmonic terms of the Fourier series. Because each Fourier series has 2*N* + 1 parameters, *n* represents the number of 1 ... 2*N* + 1.

The constraints on the RS010N robot are shown in Table 3.


**Table 3.** RS010N robot joint constraints.

When the robot moved, the observation matrix *w* was obtained by recording the angle, velocity, and acceleration of the joints of the robot in *Φb*. *n* denotes the *n*th sample.

$$w = \begin{bmatrix} \Phi\_b(q\_1, \dot{q}\_1, \ddot{q}\_1) \\ \dots \\ \Phi\_b(q\_n, \dot{q}\_{n'}\ddot{q}\_n) \end{bmatrix} \tag{8}$$

Adaptation function of the incentive trajectory.

$$y = \mathbb{C}ond(w) + P\tag{9}$$

where *P* denotes a penalty function and *Cond* denotes the number of conditions of the acquisition matrix. When the trajectory of the industrial robot satisfies the above constraint, *P* = 0, and *P* = 10<sup>8</sup> otherwise. The excitation of the first three joints of the RS010N was obtained by the genetic algorithm toolbox in MATLAB as shown in Figure 3.

**Figure 3.** The excitation trajectory of the first three joints of the RS010N robot. (**a**) First joint motivation. (**b**) Second joint motivation. (**c**) Third joint motivation.

#### **4. SEDSABSO for the Identification of Industrial Robot Dynamics Parameters**

#### *4.1. Improving Global Optimal Search*

Owing to the large number of dynamic parameters of the industrial robot that needed to be identified, their solution was a high-dimensional function optimization problem with certain constraints. The BSO used to obtain the kinetic parameter identification of industrial robots easily falls into the premature phenomenon, so the global optimal solution search mode of BSO and the iterative change form of the number *N* of tennies in the iterative process of BSO are improved in this paper.

#### BSO

BSO was proposed by Wang et al. based on the PSO and the Aspen whisker algorithm. The position of a body in BSO represents a feasible solution to the optimization problem. The formula to update BSO is:

$$
\lambda \mathbf{x}\_{i\rangle}(t+1) = \mathbf{x}\_{i\rangle}(t) + \lambda \mathbf{v}\_{i\rangle}(t+1) + (1-\lambda)\mathbf{\tilde{g}}\_{i\rangle}(t) \tag{10}
$$

where *t* denotes the number of iterations of the skink, *i* denotes the number of skinks in the skink population, *j* denotes the dimensionality of the problem to be solved, *λ* denotes *xij*(*t* + 1) denotes the position of skink *i* in the skink population at *t +* 1 iterations, ξ*ij*(*t*) denotes the position of the *i*th skink moving autonomously, and *vij*(*t* + 1) denotes the velocity of the *i*th skink. It is expressed as:

$$w\_{i\bar{j}}(t+1) = wv\_{i\bar{j}}(t) + c\_1r\_1\left(p\_{i\bar{j}}(t) - \mathbf{x}\_{i\bar{j}}(t)\right) + c\_2r\_2\left(g\_{i\bar{j}}(t) - \mathbf{x}\_{i\bar{j}}(t)\right) \tag{11}$$

*w* denotes inertial weight, *pij*(*t*) denotes the optimal position of individual *i* after *t* iterations of the Aspen population, *xij*(*t*) denotes the given position of *i* after *t* iterations, *gij*(*t*) denotes the optimal position of the Aspen population after *t* iterations, *c*<sup>1</sup> and *c*<sup>2</sup> denote adjustment factors, and *r*<sup>1</sup> and *r*<sup>2</sup> denote random numbers within the interval [0,1]

$$\xi\_{ij}(t+1) = \delta(t) \* V\_{ij}(t) \* \operatorname{sign}\left(f\left(\mathbf{x}\_{lj}(t)\right) - f\left(\mathbf{x}\_{rj}(t)\right)\right) \tag{1}$$

In Equation (12), *sign*(.) denotes the sign function, *ξij*(*t* + 1) denotes the *t* + 1th motion position searched by Aspen autonomously, *δ*(*t*) denotes the step size of the Aspen individual after *t* iterations, and *f*(*xlj*(*t*)) and *f*(*xrj*(*t*)) denote the fitness values of the left and right whiskers of the Aspen, respectively.

The positions of the left and right whiskers of the Aspen are indicated by:

$$\begin{array}{l} \varkappa\_{l\bar{j}}(t+1) = \varkappa\_{l\bar{j}}(t) + V\_{i\bar{j}}(t)d\_0/2\\ \varkappa\_{r\bar{j}}(t+1) = \varkappa\_{r\bar{j}}(t) - V\_{i\bar{j}}(t)d\_0/2 \end{array} \tag{13}$$

*xlj*(*t* + 1) denotes the position of the left whisker in iteration *t* + 1, *xrj*(*t* + 1) is the *t* + 1 right whisker position, and *d0* denotes the distance between the whiskers.

#### *4.2. Improving BSO*

4.2.1. Improved Global Optimal Solution

The global optimal particle *Pg* of BSO is randomly perturbed. Then, the new solution after perturbation is selected by the metropolis criterion of SA to improve the global convergence of BSO.

(1) Simulated annealing algorithm

The simulated annealing algorithm controls the probability of jumping out of the local optimum by setting the temperature. The algorithm evaluates whether to jump out of the worse solution by the metropolis criterion. The procedure is as follows:

a. The algorithm is initialized, and an initial solution *x* is randomly generated as the optimal solution.


Stochastic behavior is a search behavior that can improve the diversity of algorithmic populations and is widely used in a variety of intelligent algorithms, such as in the foraging behavior within the fish swarm algorithm to search for food and the wandering behavior of the wolf pack algorithm to sense the scent of prey in the air in a random pattern [17,18]. Inspired by this, we incorporated stochastic behavior into BSO, and the optimal particles in it were randomly perturbed during each iteration to enhance its ability to escape from the local optimal solution to improve global convergence:

$$P\_{\mathcal{S}}^{\
m{\rm env}} = P\_{\mathcal{S}} \ast \text{ (1} + \text{rand)} \tag{14}$$

where *Pg new* denotes the new solution after the random perturbation of the optimal particle of the Aspen swarm, and *rand* represents a random number from zero to one. The fitness value of the optimal particle *pg* of BSO was compared with that of the global optimal particle *Pg new* for random perturbation; then, the new solution was accepted according to the metropolis criterion of SA. BSO performed *q* random perturbations of its global optimal solution in each iteration.

#### 4.2.2. Improved Iterative Approach to N Particles of the Aspen Swarm

Each skyline of BSO represents a potential solution to the problem. Assuming that the number of skylines in each iteration of BSO is *N*, the maximum number of its iterations is *M*axdt, and the duration of the iteration of a single skyline is *t*, the total complexity of iterations of BSO is *N*∗*M*axdt∗*t* [19]. In solving a high-dimensional function problem, similar to determining the dynamic parameters of the industrial robot, the number of particles N at the moment of each iteration of BSO is constant. BSO is based on PSO, and a larger number of particles in PSO yields a higher accuracy of the convergence, but the time needed for convergence also increases.

For example, in previous work, we proposed a two-stage dynamic PSO that changed the number of particles in the swarm by linearly reducing the number of iterations. The experimental results showed that the accuracy of convergence of the algorithm was similar to that of the classical PSO algorithm, which also reflects a side-effect of the latter. The algorithm proposed here also had this property of a high initial accuracy of convergence [20]. Inspired by this, the number *N* of skinks in the skink swarm was changed dynamically. The number was larger in the early stage of the skink swarm algorithm and smaller in the later stage, such that the total number of skinks did not differ by much and the gap in the durations of their iterations was not large. We fully exploited the high efficiency of the early iterations of the skink swarm algorithm for the accuracy of its solution. The exponentially decreasing rate of curtailment was used for the number of cows *N* of the Aspen swarm algorithm by drawing on the idea proposed by Wang et al. [21]. That is, the particle swarm algorithm used a larger number of particles to search in the initial stage. With the increase of the number of iterations, the exponentially decreasing number of particle was used to reduce the Aspen number, so as to improve the search efficiency, and the improved iterative formula for the Aspen number is as follows:

$$N = round(e^a \* N\_{start})\tag{15}$$

$$\alpha = t \frac{\ln f\_{\text{max}} - \ln f\_{\text{min}}}{\text{Maxdt}} - \ln f\_{\text{max}} \tag{16}$$

where *Nstart* denotes the initial number of skinks of the improved skink herd, *i* denotes the number of iterations of the skink herd algorithm, *Maxdt* denotes the maximum number of iterations, and *fmax* and *fmin* are the respective maximum and minimum values of the search factor set used to control the number of skinks *N*.

#### *4.3. Robot Dynamic Parameter Identification Based on SEDSABSO*

The SEDSABSO process is illustrated below:

From Figure 4, the flow of SEDSABSO for industrial robot dynamics parameter identification is:


**Figure 4.** SEDSABSO algorithm flow chart.

#### **5. Simulation Experiments and Results Analysis**

#### *5.1. Adaptation Function*

The sum of the absolute values of the difference between the torque of joint recognition and the theoretical torque of the RS010N robot when sampled *i* times was used as the objective function, as in Equation (17)

$$\sum\_{i=1}^{N} abs(f\_k \mathbf{1}(i) - f \mathbf{1}(i)) + abs(f\_k \mathbf{2}(i) - f \mathbf{2}(i)) + abs(f\_k \mathbf{3}(i) - f \mathbf{3}(i)) \tag{17}$$

*i* = 1, 2, ... , *N* denotes the *ith* sample of the robot, *fk*1(*i*), *fk*2(*i*), and *fk*3(*i*) denote values of the first, second, and third moments of the joints values for the *i* robot recognition *f* 1(*i*), *f* 2(*i*), *f* 3(*i*) denotes the theoretical value of the three joint moments sampled by the robot.

#### *5.2. Analysis of Experimental Results*

#### Experimental Parameter Setting

To improve the accuracy of identification of the kinetic parameters, the individual algorithms were initialized uniformly by using a small interval [22]. Then, the dynamic parameters of the robot were solved for by using the SEDSABSO algorithm, BSO algorithm, and linearly decreasing particle swarm algorithm (LDWPSO), where the number of particles *N* was set to 50 for the BSO and LDWPSO algorithms. The initial population *Nstart* for SEDSABSO was set to 90, the learning factor was *c*<sup>1</sup> = *c*<sup>2</sup> = 2.0, and the inertial weight *w* was 0.8. The maximum number of iterations was set to 600. The initial temperature *T* for simulated annealing was set to 10,000, the decay scale of the annealing coefficient was set to 0.93, the step *δ*(*t*) was set to a constant value of 0.05, the adjustment factor *λ* was 0.9, and the ratio of the Aspen step to the distance between the whiskers, *c,* was set to two. After many trials, the optimal number of particle perturbations for the population was set to *q* = 3, *fmax* was 0.9, and *fmin* was 0.3. The iteration curves of the Aspen number for SEDSABSO and BSO are shown below:

Figure 5 shows that the number of particles *N* was constant during iterations of the BSO algorithm. Assuming that the maximum number of iterations was *M*axdt and the duration of iteration of each particle was *t*, the total duration of iterations of BSO was *N*∗*M*axdt∗*t*, and *N*∗*M*axdt is equal to Area1+Area2 in the figure. For the SEDSABSO algorithm, the total iteration time was approximately the product of Figure 5. For the SEDSABSO algorithm, the total iteration time was approximately the product of the area of the class trapezoid enclosed by the dash and the horizontal axis and the time *t*, i.e., (Area3 + Area1) ∗ *t*, and since the area of Area3 was not much different from that of Area2, the computational complexity of the improved Amanita group algorithm was approximately the same as that of the standard Amanita group algorithm. ѽ ѽ ѽ

**Figure 5.** SEDSABSO and BSO Variation curve of number *N* of longicorn beetles.

Figure 6 below shows the adaptation iteration curves of each algorithm on an Intel(R) Core(TM) i7-8550U main frequency 4.00 GHz computer with 600 iterations through Matlab 9.1:

**Figure 6.** Iterative curve of the algorithm.

From Figure 6, we can see that the SEDSABSO algorithm proposed in this paper converges faster and has smaller fitness values than the BSO and LDWPSO algorithms. It is clear that the SEDSABSO algorithm proposed in this paper performs best.

In order to further verify the advanced level of the SEDSABSO algorithm, this paper ran each of the above three algorithms ten times and found the average fitness and average time as follows in Table 4:

**Table 4.** Comparison of average fitness values.


From Table 4, we can see that the average fitness values of the three algorithms after ten runs were only 0.95 for SEDSABSO, 1.12 for BSO, and 1.18 for LDWPSO, further confirming the stability and effectiveness of the proposed SEDSABSO algorithm. It can be seen that the average time of SEDSABSO was better than that of BSO, further proving the analysis in Figure 5. The difference with LDWPSO was not significant, and combined with the average fitness value, it can be seen that SEDSABSO performed optimally.

The linearized kinetic parameters of the RS010N industrial robot had a total of 21 minimum parameter sets, and the results of the kinetic parameters identified by SEDSABSO, BSO, and LDWPSO are shown in Table 5 below; the absolute errors between the identified kinetic parameters and the actual parameters were also obtained.


**Table 5.** Identification of dynamic parameters.

Taking the second parameter *F*1*<sup>c</sup>* of the set of minimum parameters of the industrial robot in Table 5 as an example, the absolute errors identified by BSO and LDWPSO were −0.35 and 0.2, respectively, while the absolute error of SEDSABSO was −0.01. On the whole, the SEDSABSO algorithm yielded the smallest error in identifying the kinetic parameters of the industrial robot, while BSO and LDWPSO The kinetic parameters identified by the SEDSABSO algorithm deviated from the theoretical kinetic parameter values, which further reflects the superiority of the SEDSABSO algorithm. The maximum absolute error in the kinetic parameters as identified by the SEDSABSO algorithm was −0.39. The overall error in identification was thus small. The SEDSABSO algorithm can thus identify the kinetic parameters of the robot.

#### **6. Conclusions**


**Author Contributions:** Methodology, B.K.; Software, B.K.; Validation, B.K.; Writing—original draft, B.K.; Writing—review & editing, B.K., S.G. article chart editor, B.K. and D.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the national keyR&D plan project (2017YFB1301000).

**Data Availability Statement:** The data used in this study were self-tested and self-collected during the test. As the control method in this paper is still being further optimized, the data cannot be shared at present. Therefore, data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Maurizio Ruggiu † and Pierluigi Rea \*,†**

Department of Mechanical, Chemical and Materials Engineering, University of Cagliari, Via Marengo, 2, 09123 Cagliari, Italy; maurizio.ruggiu@unica.it


**Abstract:** This paper fits into the field of research concerning robotic systems for rehabilitation. Robotic systems are going to be increasingly used to assist fragile persons and to perform rehabilitation tasks for persons affected by motion injuries. Among the recovery therapies, the mirror therapy was shown to be effective for the functional recovery of an arm after stroke. In this paper we present a master/slave robotic device based on the mirror therapy paradigm for wrist rehabilitation. The device is designed to orient the affected wrist in real time according to the imposed motion of the healthy wrist. The paper shows the kinematic analysis of the system, the numerical simulations, an experimental mechatronic set-up, and a built 3D-printed prototype.

**Keywords:** parallel robots; mechatronics; motion simulation; mirror therapy

#### **1. Introduction**

Limb rehabilitation using robotic devices is an innovative form of rehabilitation based on interactions between the patient and the device. These systems augment the rehabilitation outcomes in neurologic disorders, such as stroke and multiple sclerosis. Robotic therapy provides high-intensity and repeated training and it can be used as an effective complement to standard rehabilitation from the beginning of a therapy [1,2]. The results show in all cases that patients who received the robotic therapy in addition to conventional therapy have greater reductions in motor impairment [3–6]. The robotic rehabilitation systems can be classified according to the design, depending on whether it operates as end-effector or as an exoskeleton. Another classification is between proximal and distal robots. Proximal robots are used to move the shoulder joint and the elbow joint; distal robots are rehabilitation robots for fine motor skills. They are used to train the hand and the fingers. Finally, rehabilitation robots can be classified as unilateral and bilateral. Unilateral devices use only the paralyzed limb for rehabilitation tasks, whereas bilateral robots use both limbs, the paralyzed and the healthy one [7].

There are numerous successful implementations of robotic rehabilitation devices. One of the first examples was the MIT-MANUS, which had a major impact on neurorehabilitation. The MIT-MANUS can move and guide the upper limb and record the trajectory, the velocity and the force of movement [8], indeed, this application can be classified as proximal and unilateral. A similar concept was described in [9]. An approach for robotic rehabilitation was used with MIME (Mirror Image Movement Enabler) [10], using an industrial robot to train the arm. A pneumatic system is the RUPERT (Robotic Upper Extremity Repetitive Therapy), which is an exoskeletal type. It consists of four pneumatic muscles, [11]. MACARM (Multi-Axis Cartesian-based Arm Rehabilitation Machine) is a planar cable-driven system designed for the rehabilitation of the upper limb [12]. The use of this kind of technology was further used in [13,14], although for correct operations, several issues for the design and modeling should be considered [15–17]. NeReBot and then MariBot are cable-driven systems with proximal, unilateral end-effector devices [18,19]. A cable-driven robot was proposed for mirror therapy in [20] and a real-

**Citation:** Ruggiu, M.; Rea, P. Development of a Mechatronic System for the Mirror Therapy. *Actuators* **2022**, *11*, 14. https://doi.org/10.3390/ act11010014

Academic Editors: Marco Carricato and Edoardo Idà

Received: 1 December 2021 Accepted: 3 January 2022 Published: 5 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

time two-axis mirror robot system was developed for functional recovery of hemiplegic arms in [21].

In this paper, we present the design, analysis and experimental set-up of a system based on the mirror therapy concept suited to wrist rehabilitation. Mirror therapy is an effective occupational therapy for functional recovery of a hemiplegic arm after stroke [22–24]. It can facilitate brain neuroplasticity through activation of the sensorimotor cortex. The standard approach uses a simple mirror and the individual sits orthogonal to the mirror. The affected limb is positioned behind the mirror, so that it blocks the view and shows the non-paralyzed limb. Watching the mirror, the motion of the master limb is ideally projected to the paralyzed limb (the slave). The mirroring creates the illusion where it looks like the paralyzed limb would do the same movement as the non-paralyzed one. With this visual illusion, damaged nerve connections in the brain should be stimulated to make reconnections [25]. Moreover, mirror therapy has been used for chronic pain [26,27], and, as it was introduced by Vilayanur S. Ramachandran [28], as a therapy against phantom pain.

In our implementation, the system is designed to orient the affected wrist (slave) in real time according to the imposed motion of the healthy wrist (master).

The rest of the paper is organized as follows: In Section 2 a description of the system is provided, Section 3 presents the kinematic analysis of the mechanisms, Section 4 shows the numerical results obtained from the simulations whilst in Section 5 we describe the experimental set-up. Finally, in Section 6 the conclusions are drawn.

#### **2. Description of the System**

The system is basically composed of two units, which are called master and slave, the master is interacting with the full functional upper limb, while the slave is devoted to the affected upper limb. The two mechanisms are different because they must have different functions. The master is designed and built as a three-axis gimbal. Figure 1 illustrates a view of the master. In a typical gimbal, there are two or three motors on the system with the aim to prevent or eliminate vibration or locate an end-effector in space [29,30]. The basic aim of a gimbal is to minimize the vibration in video recording devices, and creating a reverse motion in the opposite direction of the vibration. The reverse motion can be provided by using an inertial measuring unit (IMU) sensor, which is placed on the camera and detects the camera movements and reports the motion to three servomotors positioned in line with the camera lens. The IMU detects the relative pose of the camera according to the ground, and based on the predetermined optimum position, the deviation between the two is evaluated. Then, an electronic board receives and processes data from the IMU and then transmits the information to the servomotors of the gimbal, which provides smooth motion. Thus, the servomotors that produce the opposite movement of the camera allows obtaining a smooth image. We have used the same concept to design a gimbal for tracking the orientation of the full-functional upper limb. Instead of using the IMU, we have chosen encoders. The motivation is that the orientation of the hand must be tracked and measured, instead of its angular velocity and acceleration. Nevertheless, in future applications, we may consider the use of an IMU.

The slave is a spherical parallel mechanism [31,32]. Figure 2 illustrates a view of the slave. It consists of a fixed base and an end-effector, hereafter referred to as the joystick, connected to each other by three identical limbs, each with an RRUR kinematic chain (R stands for revolving joint, U stands for universal joint formed by two concurrent R joints, P stands for prismatic joint and S stands for spherical joint) [33]. In each limb the first three *R* joints have parallel axes forming a planar chain whilst the last two *R* joints are perpendicular to each other, intersecting in the center of rotation of the joystick. Only the *R* pairs connected to the frame are actuated. The motors torques allow the joystick to fully rotate about the center of rotation. The mechanism is decoupled and not-overconstrained.

In literature, there are numerous applications of the parallel mechanisms for rehabilitation tasks. Some ankle rehabilitation devices proposed have the 3-RRR [34], 2- RRR/UPRR [35], 2-UPS/RRR [36], 3-RRS [37] geometries and to fit more closely the ankle

motion Zhang et al. exploited a more complex parallel geometry [38]. There are, also, parallel geometries exploited for upper limb, wrist rehabilitation devices as those proposed in [39,40].

In summary, the choice of the three-axis gimbal geometry is motivated by its simplicity. In fact, the simple kinematics allows us to obtain the joystick orientation straight from the encoders measurements. Additionally, there are no strict requirements in terms of dynamics as the motion is imposed by the healthy hand. Conversely, the slave mechanism must reduce the inertial effects as the affected hand follows the driven motion and it has to be accurate and repeatable in posing the joystick. For these reasons, the parallel architecture appeared to be appropriate as it presents light moving links with the motors fixed at the base. Furthermore, the closed chain geometry may allow us to have high accuracy and·repeatability.

**Figure 1.** The master mechanism.

**Figure 2.** The slave mechanism.

#### **3. Kinematics**

The master mechanism is moved by the patient healthy arm with the rotations measured by encoders. A forward position kinematics allows us to obtain the orientation of the joystick. Because of the numerical efficiency, a quaternion parametrization was used such that the rotation of the joystick can be expressed as the compositional rotation of the intrinsic Tait–Bryan angles, *ψ*, *θ*, *φ*:

$$q = q\_3 \odot q\_2 \odot q\_{1\prime} \tag{1}$$

with

$$\begin{aligned} q\_1 &= \cos(\frac{\psi}{2}) + \sin(\frac{\psi}{2}) \mathbf{k}, & yaw \\ q\_2 &= \cos(\frac{\theta}{2}) + \sin(\frac{\theta}{2}) \mathbf{j}, & pitch \\ q\_3 &= \cos(\frac{\phi}{2}) + \sin(\frac{\phi}{2}) \mathbf{i}. & roll \end{aligned}$$

Equation (1) indicates a rotation *q*<sup>1</sup> followed by rotation *q*<sup>2</sup> and followed by rotation *q*<sup>3</sup> as shown in Figure 3.

**Figure 3.** Tait–Bryan angles in the master mechanism.

The slave mechanism drives the patient affected arm by three servo-motors. An inverse kinematics is needed to obtain the motors rotations from the orientation of the joystick *q*.

With reference to the Figure 4, the equations required to solve the inverse position kinematics are obtained considering the *i*th. of the 3-RRUR mechanism:

$$\mathbf{a}\_i + \mathbf{b}\_i + \mathbf{d}\_i = \mathbf{s}\_{i\prime} \tag{2}$$

$$\mathbf{v}\_i^T \mathbf{w}\_i = 0.\tag{3}$$

The solving procedure is: (*a*) to obtain the orientation *λ<sup>i</sup>* of **d***<sup>i</sup>* = *d***w***<sup>i</sup>* by solving Equation (3), (*b*) to obtain the motors rotations by solving Equation (2).

**Figure 4.** Slave mechanism at home pose.

Under the rotation *q* = *q*<sup>0</sup> + *q*1**i** + *q*2**j** + *q*3**k** from the joystick of the master mechanism, **w***<sup>i</sup>* are obtained as (Figure 5):

$$\begin{aligned} \mathbf{w}\_1 &= e\_1 \mathbf{j} e\_{1'}^\* \\ \mathbf{w}\_2 &= e\_2 \mathbf{k} e\_{2'}^\* \\ \mathbf{w}\_3 &= e\_3 \mathbf{i} e\_{3'}^\* \end{aligned}$$

where *e*∗ *<sup>i</sup>* is the *i*th. conjugate quaternion and

$$\begin{aligned} e\_1 &= e\_{01} + e\_{11} \mathbf{i}\_{\prime} \\ e\_2 &= e\_{02} + e\_{22} \mathbf{j}\_{\prime} \\ e\_3 &= e\_{03} + e\_{33} \mathbf{k}\_{\prime} \end{aligned}$$

*e*1,*e*2,*e*<sup>3</sup> are the Euler parameter quaternions that represent the unknown rotations *λ<sup>i</sup>* of **w***<sup>i</sup>* around the first three axes of each limb. On the other hand the axes **v***<sup>i</sup>* connected to the joystick are rotated by *q*.

$$\begin{aligned} \mathbf{v}\_1 &= q \mathbf{k} q^\* \\ \mathbf{v}\_2 &= q \mathbf{i} q^\* \\ \mathbf{v}\_3 &= q \mathbf{j} q^\* \end{aligned}$$

Equation (3) leads to:

$$\begin{aligned} &(2e\_{11}^2 - 1)(2q\_0q\_1 - 2q\_2q\_3) - 2e\_{01}e\_{11}(2q\_1^2 + 2q\_2^2 - 1) = 0, \\ &(2e\_{22}^2 - 1)(2q\_0q\_2 - 2q\_1q\_3) - 2e\_{02}e\_{22}(2q\_2^2 + 2q\_3^2 - 1) = 0, \\ &(2e\_{33}^2 - 1)(2q\_0q\_3 - 2q\_1q\_2) - 2e\_{03}e\_{33}(2q\_1^2 + 2q\_3^2 - 1) = 0. \end{aligned} \tag{4}$$

Equation (4) with the Euler parameters normalized equations, namely *e*0*<sup>i</sup>* <sup>2</sup> + *eii* <sup>2</sup> = 1, form two equations with two unknowns for each limb that can be easily solved. The sought angles *λ<sup>i</sup>* are then obtained as: *λ<sup>i</sup>* = 2*atan*(*eii*/*e*0*i*). It is worth noting that geometrically the solutions are obtained as intersection points of a conics and a circle in the Euler parameters plane {*e*0*i*,*eii*}. An example is shown in Figure 6. It may be noted that, because of the

symmetry of the solutions, a couple of points provides the same angle leading to only two distinct solutions.

**Figure 5.** Slave mechanism in an arbitrary pose.

**Figure 6.** The conics and the circumference in the Euler parameters plane {*e*01,*e*11}.

The rest of the solution procedure consists of solving the inverse kinematics of a three-link planar manipulator, namely Equation (2), with *λ<sup>i</sup>* = 3 ∑ *i*=1 *θ*1*<sup>i</sup>* as known values:

$$ac\_{\theta\_{1i}} + bc\_{(\theta\_{1i} + \theta\_{2i})} = p\_{y\_i} - dc\_{\lambda\_{i'}} \tag{5}$$

$$as\_{\theta\_{1i}} + bs\_{(\theta\_{1i} + \theta\_{2i})} = p\_{z\_i} - ds\_{\lambda\_{i'}}$$

with (*pyi* , *pzi* ) coordinates of the center of rotation expressed in the reference system of the limb. The procedure is well-known [41] and for this reason it is not reported here in the detail. First, *θ*2*<sup>i</sup>* is obtained by squaring and summing up Equation (5) to have the multiple solutions corresponding to the *elbow-up* or the *elbow-down*, then *sθ*1*<sup>i</sup>* and *cθ*1*<sup>i</sup>* can be obtained as the solution of a linear system.

#### *Range of Motion*

In the light of the application under study, it is required to know the range of motion of each input link of the slave mechanism.

In doing that the planar three-link model of the *i*th. limb is considered. According to Figure 7, the maximum Euler angle *β* and the limiting values of the input angle *θ*1*<sup>i</sup>* can be obtained geometrically. Firstly, the limit of the counterclockwise input rotation *θ*1*iu* is considered. *θ*1*iu* is reached when the first two links of the limb are aligned. This configuration represents an inverse singular configuration for the mechanism (serial singularity) that is avoided in practice:

$$\begin{aligned} (a+b)c\_{\theta\_{1i\_{\rm li}}} - ds\_{\beta} &= p\_{y\_{i'}} \\ (a+b)c\_{\theta\_{1i\_{\rm I}}} + dc\_{\beta} &= p\_{z\_{i'}} \end{aligned} \tag{6}$$

Solution of Equation (6) is straightforward. By squaring and summing up the equations an equation of the form *As<sup>β</sup>* + *Bc<sup>β</sup>* + *C* = 0 is obtained and solved by the half-tan method. Eventually, we obtain *θ*1*iu* = *atan*2(*pzi* + *dsβ*, *pyi* − *dcβ*). Once *β* is known, the limit of the clockwise input rotation *θ*1*id* is obtained by solving the Equation (7):

**Figure 7.** Geometric model for *θ*1*iu* , *θ*1*id* calculation.

A dimensional graphical synthesis was performed in order to avoid singularity configuration during the system's operations. In addition, the Jacobian matrix, which refers

to check the inverse singularity is a diagonal matrix whose condition number has been verified [33]. In particular, it is worth noting that the maximum value of the Euler angle attained by the joystick in the simulation of the system was *β*˜ = 35◦. For this value the condition number of the Jacobian matrix is about 1.7, which is still far away from singular configurations.

#### **4. Simulation of the System**

The main goal of the numerical simulations was to produce the motion mapping from the master to the slave mechanism in accordance with the results from the kinematic analysis. Further, the model allowed us to choose the motors to be used in the prototype.

The geometrical dimensions of the slave in the model are: *ai* = 72 mm, *bi* = 150 mm, *di* = 78 mm.

The simulation was carried out in order to have the same behavior (rotational motion) of the two joysticks. In particular, the motion laws of the three motors of the slave must be related to the angular configuration of the master. It is worth noting that, since the two mechanisms are different, the laws of motion will be different. The outcome of the simulation of the system has to be reproducing the same motion of the joysticks. Figure 8 shows the rotation about *z*-axis for the first 10 s and an arbitrary spherical motion for the last 10 s. Figure 9 shows the mechanisms configuration at *t* = 7.5 s when the maximum discrepancy was found between the input and the output curves while performing the basic rotation. The mapping mismatch is shown in the plot, too. For the sake of clarity, the curve is shifted by an opportune offset.

**Figure 8.** Master to slave *Rz* mapping.

**Figure 9.** Mechanism configuration at *t* = 7.5 s.

Similarly, Figures 10 and 11 show the rotation about *x*-axis in the time interval 10–20 s and the rotation about *y*-axis in the time interval 20–30 s. The corresponding mapping mismatch curves are shown as well. The mechanisms configurations at *t* = 17.5 s and *t* = 27.5 s are shown in Figures 12 and 13.

**Figure 10.** Master to slave *Rx* mapping.

**Figure 11.** Master to slave *Ry* mapping.

**Figure 12.** Mechanism configuration at *t* = 17.5 s.

It may be noted from the simulation plots that, because of the decoupled nature of the slave, the map is close to an identity while doing the basic rotations. A slight mismatch, i.e., 4.5◦, was found when the motor rotation approaches its maximum *<sup>θ</sup>Mi* <sup>→</sup> ˜ *θ*1*iu* . On the other hand, because of the complex geometry and kinematics of the slave, the motor rotation laws are different from the master counterparts while performing a spherical motion. The limiting values of the motor rotations obtained from the simulation are ˜ *<sup>θ</sup>*1*iu* = ˜ *θ*1*id* = 42◦ and the maximum value of the joystick Euler angle of *β*˜ = 35◦. As expected, all values are lower than the theoretical counterparts calculated according to Equations (6) and (7), namely *θ*1*iu* = 65.7◦, *θ*1*id* = 41◦ and *β* = 42◦. Figure 14 shows the mechanisms configuration at *t* = 38.5 s when a spherical rotation is performed.

**Figure 14.** Mechanism configuration at *t* = 38.5 s.

#### **5. Experimental Set-Up of the System**

The layout of the proposed system is shown in Figure 15 following a mechatronic architecture shown in [42,43]. The input signals of the servo-system are the reference signals *θ*1, *θ*2, *θ*3-SET, which correspond to the time position laws of the gimbal. They are sent to control board, which uses a control law *Gc* for generating the signals *θM*<sup>1</sup> , *θM*<sup>2</sup> , *θM*<sup>3</sup> -SET for the motors. *Gc* implements the inverse kinematics of the slave mechanism. In the closed-loop regulator, the *θMi* -SET are compared with the feedback signals *θMi* -F/B, which are compensated by the signal conditioning *K<sup>θ</sup>* of the angular transducers connected to the motors shafts. The error  is compensated through the proportional controller *P*, such that the outputs *θM*<sup>1</sup> , *θM*<sup>2</sup> , *θM*<sup>3</sup> -REF are the command signals of the three small servomotors of the slave system. The *θMi* -OUT signals, which are related to the shaft angular position of the motors, change the configuration of the slave to follow the *θi*-SET signals of the master minimizing the error.

\*F &RQWURO/DZ IURP0DVWHUWR6ODYH

**Figure 15.** Block scheme of the control.

Figure 16 shows the overall control scheme, in which a PC, labeled (6), is used for programming and monitoring the operation and finally recording the trial. It is worth noting that after initial programming and calibration, the system is able to work and to interact automatically with the end-user. The end-user grasps the joysticks of the master and slave systems at the same time. The healthy hand plays as the master; therefore, the gimbal (1) is completely passive and follows the movements. The encoders, named as S1, S2 and S3, measure the angular configuration (orientation) of the system and their values, *θi*-SET are sent through a sensors control board (2) to the Arduino control board (3). The resulting signals *θM*<sup>1</sup> , *θM*<sup>2</sup> , *θM*<sup>3</sup> -SET are generated by the *Gc* control law and sent to the motor control board (4) to operate the motors of the slave system to move the affected hand. It is worth noting that either a mirror movement or the same of the slave can be generated, according to the programming of the system at the beginning of the trial.

The Arduino Electronic Board is used to control the small servomotors to actuate the slave, they have a power supply of 4.8 V, maximum torque 2.45 Nm and rotation speed of 2 *π*/*s*. Figure 16 shows a scheme with the basic components, they are used for the primary task of monitoring and eventually recording the orientation of the healthy and affected hands. Nevertheless, the use of force sensors on the joysticks may be considered in order to increase safety during the operation, i.e., the system stops if the grasping forces in one or both hands increase or decrease rapidly. This option will be further implemented.

The preliminary mechatronic prototype of the system is reported in Figures 17 and 18. As mentioned previously, the system is basically composed by two units, the master is interacting with the full functional upper wrist, while the slave is devoted to the affected upper wrist.

Another important outcome of the proposed system is moving towards the e-healthcare, providing a portable system, relatively at low cost, with mechatronic solutions that allow the remote monitoring and recording of the trials to be supervised even at a distance.

**Figure 16.** System control scheme.

**Figure 17.** An overview of the proposed test bed of the mechatronic system.

**Figure 18.** Top view of the proposed test bed of the mechatronic system.

#### **6. Conclusions**

We presented the design, analysis and experimental set-up of a master/slave system dedicated to wrist rehabilitation. The basic idea was to convey the motion of the healthy wrist imposed to a gimbal-like mechanism to the affected wrist moved by a mechanism with a parallel geometry. This solution combines the simplicity of the master design with the accuracy of the slave geometry. According to the analysis and simulation outcomes, the master/slave map is an identity except when the motor links rotations of the parallel mechanism approach the upper limits. The ranges of rotation either of the joystick or of the driving motors can be considered acceptable for the application proposed. The simulations allow the authors to select the driving motors and to arrange a mechatronic set-up of the entire system. Experimental tests will be the main subject of the future work.

**Author Contributions:** Conceptualization, M.R. and P.R.; software, M.R. and P.R.; formal analysis, M.R. and P.R.; investigation, M.R. and P.R.; writing—original draft preparation, M.R. and P.R.; writing—review and editing, M.R. and P.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

