*Article* **Nonlinear Intelligent Control of Two Link Robot Arm by Considering Human Voluntary Components**

**Mingcong Deng \*, Shotaro Kubota and Yuanhong Xu**

Department of Electrical and Electronic Engineering (The Graduate School of Engineering), Tokyo University of Agriculture and Technology, 2-24-16 Nakacho, Tokyo 184-8588, Japan; 1212kubota@gmail.com (S.K.); xuyuanhongdp@gmail.com (Y.X.)

**\*** Correspondence: deng@cc.tuat.ac.jp

**Abstract:** This paper proposes a nonlinear intelligent control of a two link robot arm by considering human voluntary components. In general, human arm viscoelastic properties are regulated in different manners according to various task requirements. The viscoelasticity consists of joint stiffness and viscosity. The research of the viscoelasticity can improve the development of industrial robots, rehabilitation and sports etc. So far, some results have been shown using filtered human arm viscoelasticity measurements. That is, human motor command is removed. As a result, the dynamics of human voluntary component during movements is omitted. In this paper, based on the feedforward characteristics of human multi joint arm, a model is obtained by considering human voluntary components using a support vector regression technique. By employing the learned model, a nonlinear intelligent control of two link robot arm is proposed. Experimental results confirm the effectiveness of this proposal.

**Keywords:** nonlinear intelligent control; support vector regression; feedforward control; human arm viscoelastic

#### **1. Introduction**

In recent years, in the medical and welfare fields, human resources with appropriate skills are required for treatment/surgical support for patients and long-term care for the elderly. However, the shortage of human resources due to the declining birthrate and aging population has become a problem. As one of the solutions to the above problems, it is conceivable to adopt robots as a labor force. In the future, the places where robots will be active in society will increase not only in factories, but also in facilities and general households where they have contact with humans, so it is necessary to operate robots in cooperation with humans. Therefore, it is desirable that the robot has an excellent manmachine interface, has an affinity with humans, and has the same kinetic characteristics as humans.

Multi-joint viscoelastic properties are attracting attention as an elucidation of human motion control principles. Human arm multi-joint viscoelasticity is the characteristic of the arm joint when the human arm comes into contact with the outside world. The torque that moves a human skeletal joint is generated by the difference in tension between the leading and competing muscle groups. The above tension difference is caused by the activity of muscles controlled by commands from the central nervous system. Muscle control is used not only to generate the joint torque required for exercise, but also to change the stiffness of joints during exercise and at rest. The hardness of the joints mentioned above plays an important role in stabilizing posture and interacting with the outside world [1]. For example, when a person performs a movement to move a cup, the target is mediated by the arm. In addition to interacting with objects, it is also affected by the multi-joint viscoelasticity of this arm. In other words, it is thought that humans perform the desired movement by adjusting the multi-joint viscoelasticity so that the operating environment interacts with the arm and the object. Therefore, learning the work is not only learning

**Citation:** Deng, M.; Kubota, S.; Xu, Y. Nonlinear Intelligent Control of Two Link Robot Arm by Considering Human Voluntary Components. *Sensors* **2022**, *22*, 1424. https:// doi.org/10.3390/s22041424

Academic Editors: Andrey V. Savkin and Shafiqul Islam

Received: 6 January 2022 Accepted: 9 February 2022 Published: 12 February 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/).

the work procedure, but also learning how to control the viscoelastic properties of the above-mentioned musculoskeletal system and contact with the outside world. Exploring the mechanism of adjusting the joint mechanical impedance of the musculoskeletal system is to elucidate the motor control principle of the brain that controls complex multi-joint movements, quantitatively understands the deterioration of movements caused by nerve and muscle disorders, and has human-friendly mechanical interfaces. It can be said that it is an important issue in the development of [2].

Based on the above, many studies on human arm multi-joint viscoelasticity have been conducted [3–11]. In 1998, Gomi et al. proposed a method for estimating the viscoelasticity of the human arm during exercise using a Kalman filter [12]. Furthermore, Deng et al. focused on the numerical instability caused by the Kalman filter's digit loss, and he proposed the adaptation of the UD decomposition method as a solution. Based on these previous studies, Wang proposed motion control of a robot arm considering human arm multi-joint viscoelasticity, and its effectiveness was confirmed by simulation [13]. On the other hand, there is no example of applying human arm multi-joint viscoelasticity to the control of an actual robot arm. The problem in applying it to robot arm motion control is the reproduction of voluntary motion components. The voluntary movement component is a feedforward (FF) component output from a model based on experience in the brain when exercising. In the research to actually estimate the multi-joint viscoelasticity of the human arm, the estimation is performed after removing the above voluntary movement components with a filter. Therefore, in order to apply the multi-joint viscoelasticity of the human arm to the robot arm, it is necessary to reproduce the voluntary movement component with a feedforward controller.

In this study, we focus on the reproduction of the above voluntary movement components. Specifically, we aim to design a feedforward controller that has high control performance in the control of the robot arm. Since the feedforward controller proposed in the previous research is designed based on a mechanical model, there is a concern that the control performance will deteriorate due to modeling errors when conducting actual machine experiments. Therefore, in this research, we use Support Vector Regression (SVR) [14–18], which is a kind of machine learning methods, to design the feedforward controller to reduce the modeling error and improve the tracking performance. The control system is designed based on operator theory [19–24] to compensate for the interference and uncertainty inside the controlled object that exist when controlling the robot arm. Finally, in order to confirm the effectiveness of the proposed control system, we conduct an actual machine experiment and verify its effectiveness.

In summary, the contributions of this paper are as follows: the viscoelastic properties of the multi-joint arm are measured and analyzed through experiments. Based on the characteristic of multi-joint arm viscoelastic, a controller to simulate the human body is designed, and support vector regression is used for feedforward control.

In what follows, in Section 2, as a mathematical preparation to avoid complicating this paper, Lagrange's equation of motion used for modeling and SVR theorems used in the design of feedforward control are explained. In Section 3, as a preparation for setting the problem, we introduce the human arm multi-joint viscoelasticity, robot arm modeling and the configuration of the experimental equipment used, and then raise the problem. Section 4 explains the proposed control system, where we explain the design method for the feedback controller using multi-joint viscoelasticity of the human arm, the feedforward controller based on SVR, and the stabilization controller based on the operator theory. Section 5 first describes the experimental conditions and the SVR parameter determination method based on actual machine experiments. After that, experiment is conducted to confirm the effectiveness of the feedforward controller based on SVR. In the absence of a feedforward controller, the experimental results of using a feedforward controller based on a mechanical model and the experimental results of a feedforward controller based on SVR are compared. Section 6 describes the conclusions of this study.

#### **2. Mathematical Preparation**

In this section, Lagrange's motion equation, which is necessary when deriving a mechanical model of a robot arm, is explained.

#### *2.1. Lagrange's Equation of Motion*

When deriving the equation of motion of an object, Newton's equation of motion is generally used, but in the case of a complicated mechanical system, it is often difficult to derive it by Newton's equation. Lagrange's Equations are often used in the analysis of mechanical systems because they can solve the equations of motion of complex mechanical systems more efficiently than Newton's equations of motion. However, note that the derived solution does not change from Newton's equation because it is essentially based on the same physical law as Newton's equation of motion. In this section, we derive Lagrange's equation of motion. It is divided into three sections.

#### 2.1.1. Generalized Coordinates and Nonholonomic Constraints

In order to show the dynamic behavior of mass system and rigid system, it is necessary to select physical variables appropriately. In this section, we discuss the mass point system for simplicity, but the same idea is possible for systems including rigid systems. Generally, positions are expressed using plaque points in orthogonal coordinate systems, cylindrical coordinate systems and spherical coordinate systems, but here we consider coordinates that are convenient for expressing the position (arrangement) of the entire plaque point system and defined it as generalized coordinates. A set of generalized coordinates may include parts of a Cartesian or spherical coordinate system, but may also use angles, lengths, distances, and so on.

Now, considering any geometrical arrangement that a given mass system can take, a generalized coordinate system is said to be perfect when any of these arrangements can be represented by giving coordinates. Also, the set of generalized coordinates corresponds to continuous fluctuations in some of the coordinates, whether one of them is removed and all the rest are fixed, or all but some of them are fixed. It is said to be independent when a continuous change in its geometrical arrangement can remain. Taking some generalized coordinates for a very wide class of mass and rigid systems, including robot manipulators such as akrobot, which is the subject of this study. The number of independent coordinates in it is often constant despite changes in the permissible arrangement, which is then called the degree of freedom of the system.

A mass system has less degrees of freedom when it receives a geometric constraint. If the geometric constraint can be expressed analytically by generalized coordinates and an equation that depends only on time, the constraint is nonholonomic. Now, suppose choosing (*x*1, *x*2,..., *xm*) as the complete generalized coordinate system for a mass system, the coordinate system is not independent, there are *p* holonomic constraints such as:

$$\begin{cases} h\_1(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_{\mathfrak{M}}, t) = 0 \\ h\_2(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_{\mathfrak{M}}, t) = 0 \\ \dots, \dots, \dots \dots \\ \h\_p(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_{\mathfrak{M}}, t) = 0 \end{cases} \tag{1}$$

When these constraints are independent, there are *n* = *m* − *p* independent coordinates out of *m* coordinates. The mass system has *n* degrees of freedom. Therefore, suppose that a generalized coordinate system (*q*1, *q*2,..., *qn*) that is completely and independent from the beginning is selected for the mass system of *n* degrees of freedom. In addition, if part of a complete generalized coordinate system (*x*1, *x*2, ... , *xm*) is (*q*1, *q*2,..., *qn*), then the remaining *p* of the former are determined by Equation (1). Assuming that the mass system consists of *N* mass points, it is expressed that the position vector *r<sup>i</sup>* of any mass point *mi* is determined by the generalized coordinate system (*q*1, *q*2,..., *qn*).

$$r\_i = r\_i(q\_1, q\_2, \dots, q\_n, t) = r\_i(q, t) \tag{2}$$

The velocity *v<sup>i</sup>* of this mass point is

$$
\sigma\_i = \frac{d}{dt} r\_i = \sum\_{j=1}^n \frac{\partial r\_i}{\partial q\_j} \dot{q}\_j + \frac{\partial r\_i}{\partial t}. \tag{3}
$$

This time derivative ˙*q* = (*q*˙1, *q*˙2,..., *q*˙*n*) is called general acceleration.

Since the generalized position coordinate system (*q*1, *q*2,..., *qn*) is complete and independent, the set of infinitesimal variations of coordinates (*δq*1, *δq*2,..., *δqn*) is also complete and independent. Therefore, the variation of the position *r<sup>i</sup>* of the quality point *mi* is represented by the variation *δq*<sup>1</sup> of the generalized coordinates.

$$
\delta r\_i = \sum\_{j=1}^n \frac{\partial r\_i}{\partial q\_j} \delta q\_j \tag{4}
$$

Next, assuming that the force *f<sup>i</sup>* is acting on each mass point *mi*, the increment of all the work done by *f<sup>i</sup>* is calculated under the variational *δr<sup>i</sup>* of the arrangement of the mass system.

$$\sum\_{i=1}^{N} f\_i^T \delta r\_i = \sum\_{i=1}^{N} \sum\_{j=1}^{n} f\_i^T \frac{\partial r\_i}{\partial q\_j} \delta q\_j = \sum\_{i=1}^{N} \left( \sum\_{j=1}^{n} f\_i^T \frac{\partial r\_i}{\partial q\_j} \right) \delta q\_j \tag{5}$$

The *j*th on the right side represents the force component in that direction obtained from the infinitesimal variation *δqi* of one of the generalized coordinates *qj*, and this force is called the generalized force.

$$F\_{\vec{j}} = \sum\_{j=1}^{N} f\_i^T \frac{\partial r\_i}{\partial q\_j} \tag{6}$$

Using the generalization force, (5) is expressed as,

$$\sum\_{i=1}^{N} f\_i^T \delta r\_i = \sum\_{j=1}^{n} F\_j \delta q\_j. \tag{7}$$

#### 2.1.2. Hamilton's Principle

If the momentum vector of the mass point *mi* is set as *p<sup>i</sup>* , the equation of motion of the mass point system is expressed as

$$f\_i - \frac{d}{dt}p\_i = 0\tag{8}$$

Note that with the nonholonomic constraint, this equation is redundant and can be expressed for any variation *δri*.

$$\sum\_{i=1}^{N} \left( f\_i - \frac{d}{dt} p\_i \right)^T \delta r\_i = 0. \tag{9}$$

However, since *δr<sup>i</sup>* is generally not independent, (8) is (9). Therefore, we derive Hamilton's principle from (9) and derive *n* independent equations of motion equal to *n* degrees of freedom.

$$\sum\_{i=1}^{N} f\_i^T \delta r\_i \tag{10}$$

In general, the equation represents the sum of the work done by the forces acting on all mass points in the mass system, but it is divided into a part due to conservative force and a part due to non-conservative external force. That is, the potential energy corresponding to the conservative force is *V*(*q*), the generalized force is *Fj*, and (10) is expressed as the following equation.

$$\sum\_{i=1}^{N} f\_i^T \delta r\_i = -\delta V + \sum\_{j=1}^{n} F\_i \delta q\_j \tag{11}$$

The first term on the right side is the decrease in potential energy, and the second term is the work done by the external force. Substituting (11) into (9) gives the following equation.

$$-\delta V + \sum\_{j=1}^{n} F\_i \delta q\_j - \sum\_{i=1}^{N} f\_i - \frac{d p\_i^T}{dt} \delta r\_i = 0 \tag{12}$$

Here, the third term on the left side can be rewritten as

$$-\frac{d\boldsymbol{p}\_i^T}{dt}\delta\boldsymbol{r}\_i = -\sum\_{i=1}^N \frac{d}{dt}\left(\boldsymbol{p}\_i^T\delta\boldsymbol{r}\_i\right) + \sum\_{i=1}^N \boldsymbol{p}\_i^T\frac{d}{dt}\delta\boldsymbol{r}\_i.\tag{13}$$

Also, assuming that the mass fluctuation of each mass point of the target mass point system does not occur in the time interval considered, the variation of the total kinetic energy is

$$
\delta K = \sum\_{i=1}^{N} p\_i^T \frac{d}{dt} \delta r\_i. \tag{14}
$$

Substituting (14) into the right side of (14) palce yields

$$-\frac{d\boldsymbol{p}\_i^T}{dt}\delta\boldsymbol{r}\_i = -\sum\_{i=1}^N \frac{d}{dt}\left(\boldsymbol{p}\_i^T \delta\boldsymbol{r}\_i\right) + \delta\boldsymbol{K}.\tag{15}$$

Substituting (15) into (12) yields,

$$
\delta K - \delta V + \sum\_{j=1}^{n} F\_i \delta q\_j - \sum\_{i=1}^{N} \frac{d}{dt} \left( p\_i^T \delta r\_i \right) = 0. \tag{16}
$$

This equation holds for any time interval [*t*1, *t*2] we are thinking of, so that the variation of position *δ***r***i*(*t*1) = 0 and *δ***r***i*(*t*2) = 0. This is possible because the generalized coordinate system is perfect, and when (16) is integrated over the interval [*t*1, *t*2], the fourth term on the right side disappears and the following equation holds.

$$\int\_{t\_1}^{t\_2} \left( \delta(K - V) + \sum\_{j=1}^n F\_i \delta q\_j \right) dt = 0 \tag{17}$$

Equation (17) is called Hamilton's principle for a nonholonomic mass system with *n* degrees of freedom.

#### 2.1.3. Lagrange's Equation of Motion

In order to derive an independent equation of motion equal to *n* degrees of freedom from Hamilton's theorem, we introduce a physical quantity called Lagrangian as in the following equation.

$$L = K - V\tag{18}$$

where *K* is the kinetic energy and *V* is the potential energy. Since *V* is the potential energy, it is a function of only the generalized coordinate *q*˙*j*, but *K* is a function of *q*˙*j*, *qj* and time *t*. Lagradian *L* can be written as,

$$L = L(\dot{q}, q, t). \tag{19}$$

The variation is

$$
\delta L = \sum\_{j=1}^{n} \left( \frac{\partial L}{\partial \dot{q}\_j} \delta \dot{q}\_j + \frac{\partial L}{\partial q\_j} \delta q\_j \right). \tag{20}
$$

Substituting this into (17) yields,

$$\int\_{t\_1}^{t\_2} \sum\_{j=1}^{n} \left( \frac{\partial L}{\partial \dot{q}\_j} \left( \frac{d}{dt} \right) \delta \dot{q}\_j + \frac{\partial L}{\partial q\_j} \delta q\_j + F\_i \delta q\_j \right) dt = 0. \tag{21}$$

Here, if the first term on the left side is integrated by parts, it can be seen that (21) becomes,

$$\int\_{t\_1}^{t\_2} \sum\_{j=1}^{n} \left( -\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}\_j} \right) + \frac{\partial L}{\partial q\_j} + F\_i \right) \delta q\_j dt = 0. \tag{22}$$

Since (22) must hold for any variation *δqj*, the following *n* equations must hold for the time interval *t* ∈ [*t*1, *t*2]. This is the equation of motion of the mass system described in generalized coordinates **q** = (*q*1, ... , *qn*), and is called the equation of motion of Lagrange.

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}\_j}\right) - \frac{\partial L}{\partial q\_j} = F\_j \tag{23}$$

Many of the equations of motion of mass and rigid systems with nonholonomic constraints can be derived by using Lagrange's equation of motion with the following steps.


#### *2.2. Support Vector Regression*

Support vector regression is an application of a support vector machine to a regression problem [14–18]. Support vector regression is called SVR and support vector machine is called SVM. SVM is a typical method of binary classification and has a high prediction for unknown data. It has been reported that it is possible to construct a classifier (function) with measurement accuracy. SVM uses methods such as margin maximization and kernel tricks for identification hyperplane design, and SVR is an adaptation of these methods to regression problems. Therefore, SVR has features such as high generalization performance and effectiveness even for those with non-linear input/output relationships. This section describes the procedure for deriving the regression function and the kernel function used for the regression function.

#### 2.2.1. Derivation of Regression Function

This section describes the procedure for deriving the regression function based on SVR. The regression function of SVR is expressed by the following equation.

$$f(\mathbf{x}) = \boldsymbol{\omega}^T \boldsymbol{\Phi}(\mathbf{x}) + \boldsymbol{b} \tag{24}$$

Let *f*(*x*) be the regression function, *x* be the input vector, *ω* be the regression coefficient of the feature space, *φ* be the feature function of SVR, and *b* be the bias term. The regression function is determined from the training data using an SVM-based method. In order to determine the regression function, it is necessary to derive the regression coefficient *φ*

and the bias term *b*. Let (*xi*, *yi*) be the input/output training data used to determine the function of (24). Here, the slack variable is introduced as follows. is a setting parameter.

$$\begin{aligned} \xi\_i^+ &= \begin{cases} 0 & (y\_i - f(\mathbf{x}\_i) \le \epsilon) \\ y\_i - f(\mathbf{x}\_i) - \epsilon & (y\_i - f(\mathbf{x}\_i) > \epsilon) \end{cases} \\ \xi\_i^- &= \begin{cases} 0 & (y\_i - f(\mathbf{x}\_i) \ge -\epsilon) \\ -\epsilon - y\_i + f(\mathbf{x}\_i) & (y\_i - f(\mathbf{x}\_i) < -\epsilon) \end{cases} \end{aligned} \tag{25}$$

By using the slack variables *ξ*<sup>+</sup> *<sup>i</sup>* and *ξ*<sup>−</sup> *<sup>i</sup>* , SVR is formulated as follows.

$$\min\_{\omega, b, \tilde{\xi}} \left[ \frac{1}{2} \boldsymbol{\omega}^T \boldsymbol{\omega} + \mathbb{C} \sum\_{i=1}^n \left( \tilde{\xi}\_i^+ + \tilde{\xi}\_i^- \right) \right] \tag{26}$$

It is assumed that the constraint condition

$$\begin{cases} y\_i - \omega^T \phi(\mathbf{x}\_i) - b & \le \epsilon + \xi\_i^+, \quad i = 1, \dots, n \\ \omega^T \phi(\mathbf{x}\_i) + b - y\_i & \le \epsilon + \xi\_i^-, \quad i = 1, \dots, n \\ \xi\_i^+, \xi\_i^- & \ge 0, \qquad i = 1, \dots, n \end{cases} \tag{27}$$

is satisfied. Here, *C* is the setting parameter and *n* is the number of training data. (26) maximizes the following objective function by introducing the Lagrange multiplier *λ*<sup>+</sup> *<sup>i</sup>* , *λ*<sup>−</sup> *i* , *μ*+ *<sup>i</sup>* , *μ*<sup>−</sup> *i* .

$$\begin{split} L\_{\mathcal{P}} &= \frac{1}{2} \boldsymbol{\omega}^{T} \boldsymbol{\omega} + \mathbb{C} \sum\_{i=1}^{n} \left( \boldsymbol{\xi}\_{i}^{+} + \boldsymbol{\xi}\_{i}^{-} \right) \\ &- \sum\_{i=1}^{n} \left( \mu\_{i}^{+} \boldsymbol{\xi}\_{i}^{+} + \mu\_{i}^{-} \boldsymbol{\xi}\_{i}^{-} \right) \\ &- \sum\_{i=1}^{n} \lambda\_{i}^{+} \left( \boldsymbol{\epsilon} + \boldsymbol{\xi}\_{i}^{+} - \boldsymbol{y}\_{i} + \boldsymbol{\omega}^{T} \boldsymbol{\phi}(\mathbf{x}\_{i}) + b \right) \\ &- \sum\_{i=1}^{n} \lambda\_{i}^{-} \left( \boldsymbol{\epsilon} + \boldsymbol{\xi}\_{i}^{-} + \boldsymbol{y}\_{i} - \boldsymbol{\omega}^{T} \boldsymbol{\phi}(\mathbf{x}\_{i}) - b \right) \end{split} \tag{28}$$

Since the optimal solution is the point where the partial derivative of *Lp* with respect to *ω*, *b*, *ξ*<sup>+</sup> *<sup>i</sup>* , and *ξ*<sup>−</sup> *<sup>i</sup>* becomes 0, the following equation holds for the optimal solution.

$$\frac{\partial L\_p}{\partial \omega} = \omega - \sum\_{i=1}^n (\lambda\_i^+ - \lambda\_i^-) \phi(\mathbf{x}\_i) = 0 \tag{29}$$

$$\frac{\partial L\_p}{\partial b} = \sum\_{i=1}^n \left(\lambda\_i^- - \lambda\_i^+\right) = 0 \tag{30}$$

$$\frac{\partial L\_p}{\partial \xi\_i^+} = \mathbb{C} - \lambda\_i^+ - \mu\_i^+ = 0 \tag{31}$$

$$\frac{\partial L\_p}{\partial \mathcal{J}\_i^-} = \mathbb{C} - \lambda\_i^- - \mu\_i^- = 0 \tag{32}$$

from (29), *ω* is

$$
\omega = \sum\_{i=1}^{n} (\lambda\_i^+ - \lambda\_i^-) \phi(\mathbf{x}\_i) \tag{33}
$$

Therefore, (24) is rewritten as

$$f(\mathbf{x}) = \sum\_{i=1}^{n} \left(\lambda\_i^+ - \lambda\_i^-\right) \phi(\mathbf{x}\_i)^T \phi(\mathbf{x}) + b \tag{34}$$

Substituting (29), (31) and (32) into (28) results in the following dual problem.

$$\begin{split} \max\_{\lambda\_i^-, \lambda\_i^+} L\_p &= \max\_{\lambda\_i^-, \lambda\_i^+} \left[ -\frac{1}{2} \sum\_{i=1}^n \sum\_{j=1}^n \left( \lambda\_i^+ - \lambda\_i^- \right) \left( \lambda\_j^+ - \lambda\_j^- \right) \Phi^T(\mathbf{x}\_i) \Phi(\mathbf{x}\_j) \\ &+ \sum\_{i=1}^n y\_i \left( \lambda\_i^+ - \lambda\_i^- \right) - \sum\_{i=1}^n \epsilon \left( \lambda\_i^+ + \lambda\_i^- \right) \right] \end{split} \tag{35}$$

*λ*<sup>+</sup> *<sup>i</sup>* and *λ*<sup>−</sup> *<sup>i</sup>* are

$$\sum\_{i=1}^{n} \left(\lambda\_i^+ - \lambda\_i^-\right) = 0, \quad 0 \le \lambda\_i^+, \quad \lambda\_i^- \le \mathbb{C} \tag{36}$$

Using the kernel function,

$$K(\mathbf{x}\_i, \mathbf{x}\_j) = \boldsymbol{\phi}(\mathbf{x}\_i)^T \boldsymbol{\phi}(\mathbf{x}\_j) \tag{37}$$

(35) becomes

$$\begin{split} \max\_{\lambda\_i^-, \lambda\_i^+} \left[ -\frac{1}{2} \sum\_{i=1}^n \sum\_{j=1}^n \left( \lambda\_i^+ - \lambda\_i^- \right) \left( \lambda\_j^+ - \lambda\_j^- \right) \mathcal{K} \left( \mathbf{x}\_i, \mathbf{x}\_j \right) \right. \\ \left. + \sum\_{i=1}^n y\_i \left( \lambda\_i^+ - \lambda\_i^- \right) - \sum\_{i=1}^n \varepsilon \left( \lambda\_i^+ + \lambda\_i^- \right) \right] \end{split} \tag{38}$$

The regression function is obtained from (39) using the kernel function.

$$f(\mathbf{x}) = \sum\_{i=1}^{n} \left(\lambda\_i^- - \lambda\_i^+\right) \mathbb{K}(\mathbf{x}\_{i\prime}\mathbf{x}) + b \tag{39}$$

#### 2.2.2. Kernel Function

This section introduces the kernel function used for the regression function. As mentioned above, by using the kernel function, a complicated model can be realized without explicitly calculating *φ*(*x*). However, not all functions can be used as kernel functions, and it is generally necessary to satisfy Mercer's theorem. The necessary and sufficient condition for a continuous object and square-integrable function *K*(*x*, *x* ) to have the following expansion for the eigen *λ<sup>i</sup>* ≥ 0 and the eigenfunction *φ<sup>i</sup>* is an arbitrary square-integrable function.

$$K(\mathbf{x}, \mathbf{x}') = \sum\_{i=1}^{\infty} \lambda\_i \phi\_i(\mathbf{x})^T \phi\_i(\mathbf{x}') \tag{40}$$

The following conditions are satisfied for *g*.

$$\int\_{X^{\chi}\times\chi} K(\mathbf{x}, \mathbf{x}') \, \mathbf{g}(\mathbf{x}) \, \mathbf{g}\left(\mathbf{x}'\right) d\mathbf{x} d\mathbf{x}' \ge 0 \tag{41}$$

Any function that satisfies the above theorem can be used as a kernel function. In addition, there are many kernel functions that satisfy Mercer's theorem, and the model learned by changing the kernel function is completely different. Various kernel functions have been proposed according to the application, but this section introduces the basic kernel functions that are used very frequently. There are three basic kernel functions:

$$\mathbf{K}(\mathbf{x}\_i, \mathbf{x}) = \mathbf{x}\_i^T \mathbf{x} \tag{42}$$

$$K(\mathbf{x}\_i, \mathbf{x}) = \left(\mathbf{x}\_i^T \mathbf{x}\right)^d \tag{43}$$

$$K(\mathbf{x}\_i, \mathbf{x}) = \exp\left(-\frac{||\mathbf{x}\_i - \mathbf{x}||^2}{2\sigma^2}\right) \tag{44}$$

The above equations represent a linear kernel, a polynomial kernel and an RBF kernel, respectively. The parameter *σ* in the RBF function in (44) is expressed as *yi* as the output data, *N* as the total number of data, and *y* as the average value of the output data.

$$
\sigma^2 = \frac{1}{N} \sum\_{i=1}^{N} (y\_i - \bar{y})^2 \tag{45}
$$

A linear kernel is a simple kernel function derived when *φ*(**x***i*) = **x***<sup>i</sup>* , but it is often used when a simple model is desired. Both the polynomial kernel and the RBF kernel are capable of implementing non-linear models. Since the above two kernel functions can further adjust the complexity of the model by parameters, in many cases, they are adaptively determined for the data by using cross-validation methods or the like.

#### **3. Problem Setup**

In this section, the problem setup will be described after explaining the human arm multi-joint viscoelasticity required and the robot arm used in this study.

#### *3.1. Human Arm Multi-Joint Viscoelasticity*

Human arm multi-joint viscoelasticity is a characteristic that determines the "hardness" of a person's arm joint. The torque that moves a human skeletal joint is generated by the difference in tension between the leading and competing muscle groups. The above tension difference is caused by the activity of muscles controlled by commands from the central nervous system. On the other hand, muscle control is used not only to generate the joint torque required for exercise, but also to change the "hardness" of joints during exercise and at rest. When both muscle groups between the joints have high tension, the human arm joint becomes "hard". However, when the tensions of both muscle groups are small, the human arm joint becomes easy to move. The above-mentioned "hardness" of joints has an important role in interaction with the outside world in work and stabilization in posture maintenance.

The behavior of the human musculoskeletal system is often modeled as a springdamper mass system, including the inherent characteristics of individual muscles and the characteristics of the reflex system. Since the human musculoskeletal system actually has complicated characteristics, the expression method is not unified, but in general, the coefficient of the spring characteristic is the elasticity (stiffness) and damper characteristic of the musculoskeletal system. The coefficient is called viscosity. The coefficient of change in force with respect to change in acceleration is almost determined by the inertia of the musculoskeletal system, so it is called inertia (mass). These three coefficients are collectively called the mechanical impedance parameter of the musculoskeletal system. Among the mechanical impedance parameters, the stiffness is mainly caused by the elastic properties of the muscle, which changes according to the activity level of the muscle. In the "equilibrium position control hypothesis" [25,26], there are models of the musculoskeletal motor system, utilizing the servo system composed of its elastic characteristics and reflection. It is thought that motion and external force are generated by giving the equilibrium position as a motor command. It can also be said that "learning work" is not only learning the work procedure, but also learning how to control the viscoelastic properties of the above-mentioned musculoskeletal system and come into contact with the outside world. Therefore, exploring the adjustment mechanism of the joint mechanical impedance of the musculoskeletal system is to elucidate the movement control principle of the brain that controls complicated multi-joint movements. It can be said that it is an important issue for quantitative understanding of deterioration of movement caused by nerve and muscle disorders and development of human-friendly mechanical interface.

#### *3.2. Robot Arm*

In this section, we will introduce the robot arm, which is the experimental device handled in this study, and then explain the derivation of the dynamic model. Finally, we will introduce the hardware configuration of the experimental equipment.

#### 3.2.1. Experimental Device

In this research, we conduct an experiment using the two-degree-of-freedom horizontal robot arm that imitates a human arm. Figure 1 shows the horizontal multi-joint robot. It is characterized by using lightweight aluminum for each link and driving each link by a direct drive method.

**Figure 1.** Robot arm used in this research (Mechanical part).

The self-made buffer circuit used in the experimental equipment of the robot arm is shown in Figure 2. This circuit consists of two boards, the first stage has an input connector for a rotary encoder and an output connector for connecting to a PCI board. The motor controller for the Link 2 motor is also screwed to his first stage, but the power system is separate. In the second stage, the width of the input/output voltage differs between the PCI board and the motor controller, so an amplifier circuit composed of operational amplifiers is built.

**Figure 2.** Robot arm used in this research (Electric part).

3.2.2. Mechanics Modeling

The dynamic model of the two-link robot arm that is the control target is shown in Figure 3. Lagrange's equation of motion is used to derive the dynamic model. The equations of motion of Lagrangian and Lagrange are shown below.

$$L = K - V\tag{46}$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\chi}}\right) - \frac{\partial L}{\partial \chi} = F \tag{47}$$

*K* is kinetic energy, *V* is potential energy, *x* is generalized coordinates and *F* is generalized force. The mechanical model derived from Lagrange's equation of motion shown in (47) is expressed as the following equation [27]. (48) represents link 1 and (49) represents link 2.

$$m\_{11}\ddot{\theta}\_1 + m\_{12}\ddot{\theta}\_2 + \sum\_{j=1}^{\infty} \left\{ m\_{14}(j)\vec{v}\_{2j} \right\} + f\_1 + B\_1\dot{\theta}\_1 = \tau\_1 \tag{48}$$

$$m\_{12}\ddot{\theta}\_1 + m\_{22}\ddot{\theta}\_2 + \sum\_{j=1}^{\infty} \left\{ m\_{24}(j)\ddot{v}\_{2j} \right\} + f\_2 + B\_2\dot{\theta}\_2 = \tau\_2 \tag{49}$$

The parameters used in (48) and (49) are shown in the following equations and Table 1.


*m*<sup>11</sup> = *Is* + *Ime* + *Ie* + *meL*<sup>2</sup> <sup>1</sup> + \$ *L*<sup>1</sup> 0 *ρ*1*A*1*x*<sup>2</sup> <sup>1</sup>*dx*<sup>1</sup> + \$ *L*<sup>2</sup> 0 *ρ*2*A*<sup>2</sup> ⎧ ⎨ ⎩ *L*2 <sup>1</sup> + *<sup>x</sup>*<sup>2</sup> <sup>2</sup> + ∞ ∑ *j*=1 *φ*2*j*(*x*2)*v*2*<sup>j</sup>* <sup>2</sup> +2*L*1*x*<sup>2</sup> cos(*θ*2) − 2*L*<sup>1</sup> ∞ ∑ *j*=1 *φ*2*j*(*L*2)*v*2*<sup>j</sup>* cos(*θ*2) 7 *dx*<sup>2</sup> *<sup>m</sup>*<sup>12</sup> <sup>=</sup> *Ie* <sup>+</sup> *<sup>ρ</sup>*2*A*2*L*<sup>3</sup> 2 <sup>3</sup> <sup>+</sup> \$ *L*<sup>2</sup> 0 *ρ*2 3 *x*2 <sup>2</sup> + ∞ ∑ *j*=1 *φ*2*j*(*x*2)*v*2*<sup>j</sup>* <sup>2</sup> + *L*1*x*<sup>2</sup> cos(*θ*2)−*L*<sup>1</sup> cos(*θ*2) ∞ ∑ *j*=1 *φ*2*j*(*x*2)*v*2*<sup>j</sup>* 7*dx*<sup>2</sup> *m*14(*j*) = *Ie* + \$ *L*<sup>2</sup> 0 *ρ*2*A*<sup>2</sup> + *x*2*φ*2*j*(*x*2) + *L*1*φ*2*j*(*x*2) cos(*θ*2) , *dx*<sup>2</sup> *m*<sup>22</sup> = \$ *L*<sup>2</sup> 0 *ρ*2*A*<sup>2</sup> ⎧ ⎨ ⎩*x*2 <sup>2</sup> + ∞ ∑ *j*=1 *φ*2*j*(*x*2)*v*2*<sup>j</sup>* <sup>2</sup> ⎫ ⎬ ⎭ *dx*<sup>2</sup> *<sup>m</sup>*24(*j*) = \$ *<sup>L</sup>*<sup>2</sup> 0 *ρ*2*A*2*x*2*φ*2*j*(*x*2)*dx*<sup>2</sup> *f*<sup>1</sup> = 2*ρ*2*A*<sup>2</sup> ˙ *θ*<sup>1</sup> + ˙ *θ*2 ∞ ∑ *j*=0 *v*2*jv*˙2*<sup>j</sup>* + *L*<sup>1</sup> 2 ˙ *θ*<sup>1</sup> + ˙ *θ*2 ˙ *θ*<sup>2</sup> cos(*θ*2) ∞ ∑ *j*=0 *ρ*2*A*<sup>2</sup> \$ *L*<sup>2</sup> 0 *φ*2*j*(*x*2)*dx*2*v*2*<sup>j</sup>* + 2*L*<sup>1</sup> ˙ *θ*<sup>1</sup> + ˙ *θ*2 sin(*θ*2) ∞ ∑ *j*=0 *ρ*2*A*<sup>2</sup> \$ *L*<sup>2</sup> 0 *φ*2*j*(*x*2)*dx*2*v*˙2*<sup>j</sup>* − *ρ*2*A*2*L*1*L*<sup>2</sup> 2 <sup>2</sup> <sup>+</sup> *meL*1*L*<sup>2</sup> ˙ *θ*<sup>2</sup> sin(*θ*2) *f*<sup>2</sup> = 2*ρ*2*A*<sup>2</sup> ˙ *θ*<sup>1</sup> + ˙ *θ*2 ∞ ∑ *j*=0 *v*2*jv*˙2*<sup>j</sup>* + *L*<sup>1</sup> ˙ *θ*2 <sup>1</sup> cos(*θ*2) ∞ ∑ *j*=0 *ρ*2*A*<sup>2</sup> \$ *L*<sup>2</sup> 0 *φ*2*j*(*x*2)*dx*2*v*2*<sup>j</sup>* (50)

$$\begin{split} f\_{1} &= 2\rho\_{2}A\_{2}(\dot{\theta}\_{1}+\dot{\theta}\_{2})\sum\_{j=0}^{\infty}v\_{2j}v\_{2j} + L\_{1}(2\dot{\theta}\_{1}+\dot{\theta}\_{2})\dot{\theta}\_{2}\cos(\theta\_{2})\sum\_{j=0}^{\infty}\rho\_{2}A\_{2}\int\_{0}^{L\_{2}}\phi\_{2j}(\mathbf{x}\_{2})d\mathbf{x}\_{2}v\_{2j} \\ &+ 2L\_{1}(\dot{\theta}\_{1}+\dot{\theta}\_{2})\sin(\theta\_{2})\sum\_{j=0}^{\infty}\rho\_{2}A\_{2}\int\_{0}^{L\_{2}}\phi\_{2j}(\mathbf{x}\_{2})d\mathbf{x}\_{2}v\_{2j} - \left(\frac{\rho\_{2}A\_{2}L\_{1}L\_{2}^{2}}{2} + m\_{\epsilon}L\_{1}L\_{2}\right)\theta\_{2}\sin(\theta\_{2}) \\ f\_{2} &= 2\rho\_{2}A\_{2}\left(\dot{\theta}\_{1}+\dot{\theta}\_{2}\right)\sum\_{j=0}^{\infty}v\_{2j}v\_{2j} + L\_{1}\dot{\theta}\_{1}^{2}\cos(\theta\_{2})\sum\_{j=0}^{\infty}\rho\_{2}A\_{2}\int\_{0}^{L\_{2}}\phi\_{2j}(\mathbf{x}\_{2})d\mathbf{x}\_{2}v\_{2j} \\ &+ \left(\frac{\rho\_{2}A\_{2}L\_{1}L\_{2}^{2}}{2} + m\_{\epsilon}L\_{1}L\_{2}^{2} + m\_{\epsilon}L\_{1}L\_{2}\right)\theta\_{1}^{2}\sin(\theta\_{2}) \end{split}$$

#### 3.2.3. Hardware Configuration

Figure 4 shows a schematic diagram of the experimental equipment. Motor 1 is TOYO TECHINICA DM-008D25F, motor 2 is maxon RE25 series 118752, encoder 1 is NEMICON 38H-4096-2MC and encoder 2 is NEMICON 18M-1024-2MC. The control program is written in *C*#. The command value calculated by the PC is DA-converted by PCI3521, then amplified twice by the buffer circuit and input to the servo amplifier. The voltage value is converted into a current value by the servo amplifier, and the current drives the motor to operate each link. The angle of each link is measured by capturing the number of pulses obtained from the encoder into a PC using the pulse counter board PCI6204.

#### 3.2.4. Problem Setup

Human arm multi-joint viscoelasticity is,


Among them, there are many studies on the estimation of human arm multi-joint viscoelasticity and there are few examples of application to mechanical interfaces as in the third entry. The challenge in applying it to machine interfaces is the reproduction of voluntary movement components. The voluntary movement component is a feed-forward component output from a model based on experience in the brain when exercising.

**Figure 4.** Proposed control system.

Humans suppress disturbances with feedforward controls composed of the cerebellum. It is believed that the body is controlled by a feedback controller. The voluntary movement component represents the above feedforward component, and the human arm multi-joint viscoelasticity represents the feedback component. In the research to actually estimate the multi-joint viscoelasticity of the human arm, the estimation is performed after removing the above voluntary movement components with a filter. Therefore, when applying the multi-joint viscoelasticity of the human arm to the motion control of the robot arm, it is indispensable to reproduce the voluntary motion component. In the previous research, the motion control of the robot arm by the two-degree-of-freedom control system has been introduced using multi-joint viscoelasticity. A feedforward controller based on a mechanical model has the advantage that it is easy to design if the modeling of the system is completed, but it is easily affected by modeling errors. It may cause deterioration of control performance when conducting an actual machine experiment. In fact, in this study as well, when a control experiment using a feedforward controller based on a mechanical model was conducted, good control results can not be obtained due to things that are not taken into consideration during modeling, such as the influence of the dead zone of the motor driver. Therefore, in this research, we propose a feedforward controller based on SVR, which is one of machine learning methods. SVR is an application of SVM to a regression problem, and has features such as high generalization performance and effectiveness for non-linear input/output relationships. In addition, by learning the input/output relationships of the entire system including the motor driver as training data, it is possible to create a model that includes modeling errors and parts that were not considered during modeling. In this research, we design a feedforward controller using SVR and confirm its effectiveness in actual machine experiments.

#### **4. Control System Design**

This section shows the proposed control system design method. Figure 5 shows the proposed control system. We design a two-degree-of-freedom control system by simulating the control system of the human body introduced in Section 3. The multi-joint viscoelasticity measured from humans is applied to the feedback controller *C*. The feedforward controller *F* based on SVR is also designed. In addition, in order to eliminate the influence and uncertainty due to interference inside the controlled object, the control system is designed based on the operator theory.

**Figure 5.** Block diagram of the proposed control system.

#### *4.1. Controller Design Based on Multi-Joint Viscoelasticity*

The follow-up controller *C* based on multi-joint viscoelasticity is shown by the following equation.

$$\mathbf{C}(\mathbf{e}) = \mathbf{D}\dot{\mathbf{e}} + \mathbf{R}\mathbf{e} \tag{51}$$

**R** represents an elastic matrix and D represents a viscoelastic matrix. By applying this multi-joint viscoelastic matrix to a feedback controller, we aim to reproduce the same motion as humans. The elastic matrix and the viscosity matrix are 2-by-2 matrices and can be expressed as follows.

$$\begin{aligned} \mathbf{R} &= \left( \begin{array}{cc} R\_{11} & R\_{12} \\ R\_{21} & R\_{22} \end{array} \right) \\ \mathbf{D} &= \left( \begin{array}{cc} D\_{11} & D\_{12} \\ D\_{21} & D\_{22} \end{array} \right) \end{aligned} \tag{52}$$

The multi-joint human arm viscoelasticity used in this paper is shown in Figure 6. The combined movement of Link1 and Link2 produces translational movement, as shown in Figure 7.

**Figure 6.** Multi-joint human arm elasticity and viscosity.

#### *4.2. Feedforward Controller Design Based on SVR*

This section describes the design method of the feedforward controller based on SVR. Since SVR is one of the regression analysis methods, it is necessary to select training data. PD control is performed for the target trajectory used in the experiment in this study, and the input voltage applied to the motor driver of each link, the angle of each link and the angular velocity of each link at that time are measured as training data. The training data is shown in Figure 8.

SVR learning is performed based on the results shown in Figure 8. Since the angular velocity of link 2 was greatly affected by the measurement noise, the data filtered by the RC filter was trained as training data. In addition, the controller to be designed calculates an appropriate input voltage for the target angle and angular velocity. Therefore, the input/output relationship of the data to be learned by the SVR is the angle and angular velocity at the input and the voltage at the output. Note that it is the opposite of the normal input/output relationship.

**Figure 7.** Multi-joint viscoelasticity.

**Figure 8.** Training data of Link 1 and Link 2.

*4.3. Control System Design Based on Operator Theory* 4.3.1. Elimination of Uncertainty and Interference

The nominal model *P*ˆ, which eliminates uncertainty and interference from other variables, can be expressed by the following equation.

$$\mathcal{P} = \begin{cases} \begin{array}{c} \hat{m}\_{11}\ddot{\theta}\_{1} + B\_{1}\dot{\theta}\_{1} = \mathfrak{r}\_{1} \\ \hat{m}\_{22}\ddot{\theta}\_{2} + B\_{2}\dot{\theta}\_{2} = \mathfrak{r}\_{2} \end{array} \tag{53}$$

Using this nominal model, the operators *S* and *R* are designed to eliminate the effects and uncertainties caused by interference inside the controlled object [13]. By designing the controllers *S* and *R* as *SP*ˆ = *I*, *R* = *I*, the uncertainty of the model and the influence of interference inside the controlled object can be eliminated, where *I* is an identity map. The following equation holds from the nonlinear feedback system shown in Figure 5.

$$\begin{aligned} u^\*(t) - S(y)(t) + S\hat{N}\hat{D}^{-1}(u)(t) &= R(u)(t) \\ u^\*(t) &= S(y)(t) - S\hat{P}(u)(t) + R(u)(t) \\ &= \mathcal{P}^{-1}(y) \end{aligned} \tag{54}$$

At this time, *y*(*t*) = *P*ˆ(*u*)(*t*), which shows that the uncertainty of the model and the influence of interference inside the controlled object can be eliminated. Equivalent feedback loops before and after removing uncertainty and interference are shown in Figure 9.

(**a**) Control system before eliminating interference.

(**b**) Control system after eliminating interference.

**Figure 9.** Control system before and after eliminating interference .

#### 4.3.2. Guarantee of Stability

In this section, the control system is designed based on operator theory, and the stability of the proposed control system is guaranteed. Specifically, based on operator theory, we design stabilization controllers *A* and *B*−<sup>1</sup> that satisfy the Bezout equation *AN*ˆ + *BD*ˆ = *I*. The right decomposition of the nominal plant *P*ˆ from which interference and uncertainty have been removed gives the following equation.

$$\begin{aligned} \hat{N} = \begin{cases} \,^\sharp \hat{m}\_{11} \dot{y}\_1 + B\_1 \dot{y}\_1 = \omega\_1 \\ \,^\sharp \hat{m}\_{22} \dot{y}\_2 + B\_2 \dot{y}\_2 = \omega\_2 \\ \,^\sharp \theta\_1 = y\_1, \theta\_2 = y\_2 \end{cases} \end{aligned} \tag{55}$$

$$\mathcal{D}^{-1} = \begin{cases} \ \omega\_1 = \mathfrak{r}\_1 \\ \ \omega\_2 = \mathfrak{r}\_2 \end{cases} \tag{56}$$

The parameters used are shown below.

$$\begin{aligned} \hat{m}\_{11} &= I\_s + I\_{m\varepsilon} + I\_{\varepsilon} + m\_{\varepsilon}L\_1^2 + \int\_0^{L\_1} \rho\_1 A\_1 \mathbf{x}\_1^2 d\mathbf{x}\_1 \\ &+ \int\_0^{L\_2} \rho\_2 A\_2 \left\{ L\_1^2 + \mathbf{x}\_2^2 + 2L\_1 \mathbf{x}\_2 \right\} d\mathbf{x}\_2 \\ \hat{m}\_{22} &= \int\_0^{L\_2} \rho\_2 A\_2 \mathbf{x}\_2^2 d\mathbf{x}\_2 \end{aligned} \tag{57}$$

Operator *A* is designed as follows:

$$A\hat{N} = kI\tag{58}$$

From the Bezout equation, *B* can be expressed as follows.

$$B = (1 - k)\hat{D}^{-1} \tag{59}$$

*k* is a design parameter. By constructing a nonlinear feedback system as shown in Figure 9b using operators *A* and *B*, the BIBO stability of the control system can be guaranteed.

#### **5. Experiment**

In this section, in order to confirm the effectiveness of the proposed control system, we verify it with an experimental device.

#### *5.1. Experimental Conditions*

Table 2 shows the parameters of the experimental equipment. The initial angle and target angle are the same as the angles used when the multi-joint viscoelasticity estimation was performed. Also, the sampling time was set to 0.01 s.

**Table 2.** Laboratory equipment parameters.


#### *5.2. Selection of Hyperparameters of SVR*

In this section, we select the hyperparameter *c* of SVR. There are three hyperparameters in SVR, and it is known that the regression model changes depending on the parameters. In this research, we focus on *c* among hyperparameters, experiment by changing the value of *c* step by step, and select the parameter with the best result. As the content of the experiment, for the proposed control system, only the hyperparameter *c* of the feedforward controller based on SVR is changed to 0.01, 0.05, 0.1, 1, 7, 10, and 100, and the experiment is performed. After that, the error between the experimental result and the target value is derived and evaluated by RMSE (root mean squared error). RMSE is expressed by the following equation.

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (f\_i - y\_i)^2} \tag{60}$$

Since the closer the RMSE is to 0, the smaller the error is, we select the hyperparameter *c* that minimizes the RMSE. The RMSE changes with different values of *c*, shown as Table 3. The experimental results are shown in Figure 10. Looking at the result of link 2, we can see that the control result changes depending on the value of hyperparameter *c*. When *c* = 100, it can be seen from Figure 10 that a large overshoot occurs. If *c* is set to a large value, SVR is closer to the hard margin. While it is possible to create an accurate regression model that reflects most of the training data, it has the characteristic that noise during training data

measurement is easily reflected in the regression model. Therefore, when a large value such as *c* = 100 is set, a regression model that reflects the measurement noise contained in the training data is created, and it is considered that the overshoot shown in Figure 10a occurred. On the other hand, it can be confirmed that when *c* is made smaller, the above overshoot decreases and almost disappears at the time of *c* = 0.1. It can be seen that when *c* is set to 0.1 or less, a slight delay occurs at the time of 3*s* to 4*s* due to the effect of the output from the SVR controller becoming smaller. Since the value of *c* = 0.1 is also the minimum in the evaluation of RMSE, the *c* of link 2 is set to 0.1 in this study. For link 1, *c* is selected as 0.1 in the same way as link 2.


0.01 0.0250 0.0085

**Table 3.** RMSE results with different value of *c*.

**Figure 10.** The experimental results.

#### *5.3. Experimental Results*

In this section, we conduct experiments to confirm the effectiveness of the proposed control system and introduce the results. Specifically, in order to confirm the effectiveness of the proposed SVR-based feedforward controller, it is effective to compare the experimental results of the feedforward controller based on the dynamic model without the feedforward controller with the experimental results of the proposed method. In addition, the control systems other than the feedforward controller are the same, and the performance of the feedforward controller is compared. Figure 7 shows the elastic ellipsoid of the hand calculated from the multi-joint viscoelasticity used in the feedback controller. As in the previous section, RMSE is used for comparison of experimental results, and the best result is the experimental result with the smallest RMSE. Table 4 shows the RMSE results.



Figure 11 shows a comparison of the control results between the feedforward controller without the feedforward controller and the feedforward controller based on the dynamic model. In the result of link 2, it can be confirmed that the tracking performance is improved from 2.5 s to 3 s . The result of RMSE is also smaller in the feedforward controller based on the dynamic model than in the case without the feedforward controller, and the effectiveness of the feedforward controller can be confirmed. On the other hand, for Link 1, no significant improvement is seen in the control results, and the analysis results by RMSE did not change much. The cause is thought to be the error that occurred during modeling.

**Figure 11.** Comparison of FF without controller and dynamic model with FF controller.

Figure 12 shows a comparison of the control results between the feedforward controller without the feedforward controller and the feedforward controller based on SVR. It can be confirmed that the tracking performance is improved for both link 1 and link 2. The analysis result by RMSE is also smaller in the feedforward controller based on SVR than in the case without the feedforward controller, and the effectiveness of the feedforward controller can be confirmed.

**Figure 12.** Comparison of FF without controller and SVR FF with controller.

Figure 13 shows a comparison of the control results of the feedforward controller without the feedforward controller, the feedforward controller based on the dynamic model and the feedforward controller based on the SVR. Comparing the control results of the feedforward controller based on the dynamic model and the feedforward controller based on SVR, it is confirmed that the feedforward controller based on SVR follows the target value more for both link 1 and link 2. The RMSE value is also lower in the feedforward controller based on SVR, and it can be seen that the feedforward controller based on SVR has better tracking performance numerically. It is considered that this is because it was possible to create a model closer to the experimental equipment by creating a model from the training data of the actual machine experiment using SVR.

**Figure 13.** Comparison of all control results.

Figure 14 shows the target value of the hand position coordinates and the actual output. With the proposed control system, we are able to confirm the follow-up of the hand. From the above, we are able to verify the effectiveness of the feedforward controller based on the proposed SVR method in actual machine experiments.

**Figure 14.** Multi-joint arm movement trajectory.

#### **6. Conclusions**

In this study, we propose a two-degree-of-freedom control system using multi-joint viscoelasticity, and conduct a motion control experiment of a two-link robot arm. Focusing on the feedforward controller in the two-degree-of-freedom control system, we propose a feedforward controller based on SVR, which is one of machine learning methods. Finally, the effectiveness of the proposed method is verified by an actual machine experiment. The controller is designed as a multi-joint arm like one and it is based on the characteristic of the human arm multi-joint viscoelasticity. The characteristic is analyzed from the experiment.

In the future, some intelligent control as well as adaptive control methods [28–30] can be considered to further improve the current work.

**Author Contributions:** Data curation, S.K.; Investigation, S.K.; Methodology, M.D.; Software, S.K.; Supervision, M.D.; Validation, S.K.; Writing, original draft, Y.X. 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:** Data are not publicly available due to privacy considerations.

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

#### **References**


#### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Sensors* Editorial Office E-mail: sensors@mdpi.com www.mdpi.com/journal/sensors

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8149-1