**1. Introduction**

Soft robotics has been gaining importance in the robotics research field. The intrinsic compliance and adaptable properties of this hardware are pushing them into many areas. The purpose of these technologies is to overcome some of the problems found in the current robotic platforms. These include weight, cost, versatility and more importantly, safe human-to-robot interaction.

Different soft robotics technologies have emerged. These include pneumatic muscles with rigid links [1], pneumatic materials that deform according to their strain field [2], robots with fully inflatable links [3], fully inflatable robots [4], plant-based structures [5] and many other technologies [6,7]. In particular, we are interested in tendon-driven soft robots, a bio-inspired model scheme, as those in [8–10]. However, the kinematic models, unlike rigid ones, are not yet well-understood. Given the high non-linearity and physical characteristics, several assumptions and numeric simplifications are considered to actuate. Therefore, they are not as reliable and have lower versatility in comparison with their counterpart, thus limiting their impact on robotics [11]. These drawbacks are stopping soft robotics to enter fields such as industrial robotics or manufacturing. However, where

**Citation:** Quevedo, F.; Muñoz, J.; Castano Pena, J.A.; Monje, C.A. 3D Model Identification of a Soft Robotic Neck. *Mathematics* **2021**, *9*, 1652. https://doi.org/10.3390/math9141652

Academic Editors: Mikhail Posypkin, Andrey Gorshenin and Vladimir Titarev

Received: 3 June 2021 Accepted: 9 July 2021 Published: 13 July 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/).

precision is not mandatory or humans might apply proper correction to achieve the desired goal, they have found a niche [12,13], such as rehabilitation and prosthetics.

In [8], the authors developed a mathematical model for the tendon-driven robotic arm. In particular, the authors approximated the elasticity of the tendons as a mass-less spring, given that their mass is many times smaller than the other parts (motors, gears, and loads). This enables them, on one side, to neglect coupling effects over different links, and also to assume rigid motions of the particular link to be modeled. Notice that even though the motion of the tendons is aligned with the arm motion, obtaining such a model is troublesome and requires further developments to increase the accuracy of the final result. To overcome the modeling problem for control purposes, the authors in [10] used reinforcement-learning to control the tendon-driven ACT hand synergies. This robotic hand has 24 motor-driven tendons that mimic the human hand biomechanics. Therefore, dynamic interaction with the hand skeleton results in redundant motions and other nonlinear characteristics of the hardware that should be considered if a mathematical model is required. By using reinforcement learning, the authors are able to capture the desired dynamics over a set of motions that derive into a control strategy over the 24 tendons in a reduced state-space. However, as pointed out by [9], data-driven control algorithms on tendon-based robots or soft robotics have not yet been explored, and model-based control strategies are still preferred. Considering the previous statement, the authors in [9] analyzed an autonomous learning algorithm to obtain the model of a tendon-driven leg when different stiffness is used.

From the modeling perspective, [14] proposed a finite element model (FEM) for a glove with pneumatic bending actuators. The authors worked in a two-dimensional space, neglecting the dynamic energy from the model. In further research, a black-box model identification was given by [15] for a fluidic actuator. This allows the authors to introduce the shear deformation into the model, which in their previous work was not considered. The maximum parametric variations reached 36%, which was compensated by a back-stepping controller. In the present paper, black-box data-driven modeling will be used. Therefore, the overall dynamics of the given neck will be considered and modeling results will be compared with standard linear and non-linear modeling techniques, that is, Recursive Least Square and Neural Networks. This provides an overview of alternative modeling possibilities and their implications over a non-linear system as the tendondriven neck.

In [16], a geometrical model and a two-dimensional FEM model for a soft fluidic actuator were studied. The geometrical model considered a uniform bending curvature of the link, while the FEM model showed a linear trend on the link behaviour. A new approach was proposed in [17] in order to obtain a pneumatic soft-arm 3D model, based on a constant link's bending curvature and neglecting the gravity or payload effects. Another geometrical approach for modeling a soft link is presented in [18]. In this case, the approached model neglects the effects of gravity or internal elastic forces. Regarding tendon-based robots, in the hand exoskeleton given in [19,20], each finger is considered as a three-link kinematic chain and all the joints are considered to be pure revolute. Friction and cable guide deformation were neglected. As shown in the cases above, geometric modeling requires several constraints and assumptions to reduce the model complexity, which allow the designers to cope with the complex system mechanics. Nonetheless, the black-box modeling approach will include the whole system dynamics, which led to the proposed methods in this work.

Modeling of soft robotic links is of particular interest, given the coupled dynamics that arise when actuating the robots. In this sense, different approaches are presented as well. Geometrical approaches expose their limitations when the modeling space increases. Therefore, sensor-based approaches are being used, although their results are limited to the sensor's range and capabilities. In [21], textile strain sensors are embedded into the robot structure to calculate the link deflection state and position. A similar approach is presented in [22,23]. In this last one, the authors describe the implementation of a soft hand where

each finger uses an elastic joint with an embedded piezo-electric transducer to sense the deflection of the joint. In this case, the curvature is given by the sensor but is assumed to be continuous through the link. A different sensor technique for embedded deflection modeling is presented in [24], where photo-sensing is used to determine the deflection angle of the link. A 3D modeling approach is given in [25], where embedded cameras for selfobserving of the robot configuration are synchronized to obtain the final soft body (robot) shape through learning algorithms. The named references imply additional hardware but allow for rigid modeling based on sensor information. However, their effectiveness is limited to the sensor's accuracy and the validity of the hypothesis over the soft link, such as continuous link curvature. Furthermore, to obtain the mathematical model, it is still required to neglect certain dynamics, and in some cases, only the 2D motion of the soft link is considered. In this work, the soft link is aimed to move in 3D and it is not desired to add additional hardware to the system.

This work aims to characterize the 3D motion of a tendon-driven soft neck described in [26]. An improved version of the initial version was later proposed in [27]. This new design features a soft bendable spine that replaces the central spring. The replaced structure decreases the overall weight and increases the system robustness. An inclination sensor was also introduced on the top platform for extracting the current pitch and roll angles, and enables feedback control.

Although the central structure introduced an already non-linear behaviour, as seen in [27], the new material adds other non-linearity mechanics. Some previous works already tackle system identification. Firstly, in [26] for control design, the study was limited to actuators and ignored the dynamics of the link. The resulting theoretical model is outside the standard modeling methods. As a consequence, additional methodology is required to extract a simulation and control model for the platform. An initial identification exploration on 2D was presented in [28]. In that work, which this paper is a continuation of, we identified the soft link dynamics considering the actuators and the soft link. However, only the front inclination was considered, neglecting at that time the interactions that occur when all the soft neck degrees of freedom are used.

In this work, the proposed models are improved and extended to the entire robot motion range. Set membership and Recursive least-squares identification methods are used for modeling as in [28]. As the recursive least-squares method is only valid for linear plants, the non-linear behavior will not be captured. The selected methods do not need hardware modifications nor neglected dynamics. Therefore, physical effects, such as gravity, elasticity, and plasticity will be considered by the obtained 3D models when possible. These models are compared with a neural network model identification as ground comparison. As an important contribution, no modeling technique selected relies on local deformation sensors, and they do not require additional external hardware for possible neck control considered in the future.

The remaining parts of this paper are organized as follows. In Section 2, the platform to be identified is described. Sections 3 and 4 present the different methods used for identification. In Section 5 the experimental procedures are described and Section 6 shows the resulting models. Then, in Section 7, different tests are performed for validation and comparison purposes. Finally, in Section 8, the main conclusions are discussed.

#### **2. Soft Neck Description**

The mechanism that enables soft neck operation is the central soft link, which acts as a spine. It is made with bendable material and actuated with a parallel mechanism driven by cables, which produce a tilt in the upper platform. Figure 1 shows the soft neck prototype and its parts.

The neck is composed of a base, a mobile platform, a central soft link, tendons (cables), and motors, as shown in Figure 1. All parts were built using a 3D printer, including the soft link, which weighs 100 gr (excluding motors and hardware).

Combining the actions of the three actuators, any position or orientation can be reached inside the bounded space. As the robot workspace is three-dimensional, the final rotation can be defined using three Euler angles. The *Z*-axis rotation (yaw) is neglected since it cannot change due to the configuration of the link; therefore, two rotations around the *X*- and *Y*-axes are enough to fully define the robot position and orientation. System output will be defined through an *X*-axis rotation, roll(*φ*), and a *Y*-axis rotation, pitch (*θ*), as shown in Figure 2.

**Figure 1.** Soft neck platform.

**Figure 2.** Soft neck kinematics. Orientation and inclination variables [27].

According to [29,30], the soft neck is a hyper-redundant robot. Therefore, the term degrees of freedom (DOF) is not applicable in the usual sense. Nevertheless, there is a connection between the three tendon lengths and the neck's final angular position.

The three tendon actuators are located at the base, each composed with a motor, gear, encoder, and a driver, with the characteristics shown in Table 1.


**Table 1.** Platform hardware specifications.

There exists low-level control managed by the motors' drivers, which satisfies all the platform's needs. For this reason, all system data are captured as an open-loop plan, with the actuator position and velocity as input and the inclinations roll and pitch as output.

## *2.1. Geometric Simplification*

The described plant input and output variables are coupled, defining a Multiple Inputs Multiple Outputs (MIMO) system, which makes the system identification more difficult. Fortunately, there is a way to simplify the system, decoupling these variables, allowing to analyze its behavior through simpler Single-Input Single-Output (SISO) models. Using this scheme, the inverse kinematics described in [31] will not be necessary, making the model identification easier.

The angle combination produced by each actuator effect is a final rotation that can be defined or measured by two angles. A more detailed analysis of the robot's geometry will show the connection between the inputs and these final angle outputs. Since the inputs and outputs of the system are coupled in a MIMO system of three inputs and three outputs, it is desirable to rearrange the original input scheme using a linear combination of them. In this way, three new inputs will be obtained, having a direct action with respect to the outputs. Note that the aim is not to simplify the system, but to study the effects of the different inputs in the neck output variables. To simplify this operation, we will consider the system in the resting position. In this state, the single effect of the *A*1 actuator (reducing *l*1) results in a rotation aligned with the *X*-axis of the base frame (see Figure 2). This results in an output angle directly related to the length of its tendon, and therefore, the motor position. Given that motor angles can be negative, we will consider the neck's rest position (*pitch* = 0,*roll* = 0) as the initial zero value for all tendons, resulting in positive values when tendon lengths increase, and negative otherwise. Considering just the first actuator with the index number 1, a possible equation describing pitch angle *θ* in *X* is the following:

$$
\theta\_1 = f(P\_1) \tag{1}
$$

where *θ*<sup>1</sup> is the angle contribution from the first actuator to the final pitch angle (*θ*), *P*<sup>1</sup> is the actuator input position, and *f* is a nonlinear function describing the relation between both. Although it is considered that just *P*<sup>1</sup> can change the angle *θ*<sup>1</sup> as shown in Equation (1), given the neck's nonlinear nature, the other inputs may have a tiny effect on that angle, too, but they are considered too small and will be neglected in this case. The other actuators' effects (*θ*2, *θ*3) on the final angle *θ* are the following:

$$\theta\_2 = \cos(\gamma\_2) f(P\_2) \tag{2}$$

$$\theta\_{\mathfrak{Z}} = \cos(\gamma\_{\mathfrak{Z}}) f(P\_{\mathfrak{Z}}) \tag{3}$$

Given the proposed vertical robot setup, and using the same actuators, we can assume that the functions *f* are similar. Nevertheless, a projection factor needs to be considered, which depends on the actuator's relative angle (*γ*). We can generalize the previous functions in the following equation:

$$\theta\_i = \cos(\gamma\_i) f(P\_i) \tag{4}$$

Keeping in mind that we are not modeling the system but proposing an alternative input–output scheme, we can consider these input angles additive. Again, it is not a model simplification, but an input redefinition.

$$\theta = \theta\_1 + \theta\_2 + \theta\_3 = \cos(\gamma\_{11})f(P\_1) + \cos(\gamma\_{12})f(P\_2) + \cos(\gamma\_{13})f(P\_3) \tag{5}$$

According to Figure 2, the three actuators are symmetrically arranged, and the angles are *γ*<sup>11</sup> = 0 deg, *γ*<sup>12</sup> = 120 deg and *γ*<sup>13</sup> = 240 deg. Therefore, Equation (5) results in:

$$\theta = f(P\_1) - 0.5f(P\_2) - 0.5f(P\_3) = f(P\_1) - 0.5[f(P\_2) + f(P\_3)].\tag{6}$$

This result shows how both *A*<sup>2</sup> and *A*<sup>3</sup> actuator effects on the angle *pitch* are divided by two, with an opposite direction to actuator *A*1. This leads to the first result of this approach. The *pitch* angle is defined by the length difference, being positive when *P*<sup>1</sup> is larger than 0.5(*P*<sup>2</sup> + *P*3), and negative otherwise. In the case of *P*<sup>1</sup> = *P*<sup>2</sup> = *P*3, angle *θ* = 0, leading to different robot compression states depending on the tendon lengths, form zero (*P*<sup>1</sup> = *P*<sup>2</sup> = *P*<sup>3</sup> = 0) to full compression ( *P*<sup>1</sup> = *P*<sup>2</sup> = *P*<sup>3</sup> = *Pmax*). This feature could be used to change the neck stiffness, although this is not discussed in this paper, where the soft link length is considered constant.

Now, roll angle (*φ*) is defined as the rotation around the *Y*-axis. Using the previous reasoning but projecting in the *Y*-axis (using sin(*γ*)):

$$
\phi = \phi\_1 + \phi\_2 + \phi\_3 = \sin(\gamma\_1)f(P\_1) + \sin(\gamma\_2)f(P\_2) + \sin(\gamma\_3)f(P\_3) \tag{7}
$$

In the case of *γ*<sup>1</sup> = 0 deg, *γ*<sup>2</sup> = 120 deg and *γ*<sup>3</sup> = 240 deg, Equation (7) results in:

$$
\phi = 0.8666 f(P\_2) - 0.8666 f(P\_3) = 0.8666 [f(P\_2) - f(P\_3)] \tag{8}
$$

Note that in this case, the value of the *φ* angle just depends on the difference between *P*<sup>2</sup> and *P*3, and that the *A*<sup>1</sup> actuator has no effect. Again, the angle just depends on their difference, and the compression is an average function of the tendon lengths. For the case *P*<sup>1</sup> = *P*<sup>2</sup> = *P*3, angle *φ* = 0, leading to the same previous result regarding soft link compression. Additionally, note that *θ* and *φ* angles depend on the tendon length difference, and the compression (*δ*) depends on the tendon lengths' average. Based on this, we can define the new input variables *θi*, *φ<sup>i</sup>* and *δ<sup>i</sup>* as a linear combination of the motor position inputs.

Using the results from Equations (6) and (8), and considering the link compression input as the motor positions' average, the following input redefinition is proposed:

$$\theta\_i = P\_1 - 0.5(P\_2 + P\_3) \tag{9}$$

$$
\phi\_i = 0.866(P\_2 - P\_3) \tag{10}
$$

$$\delta\_i = \frac{P\_1 + P\_2 + P\_3}{3} \tag{11}$$

Using this input redefinition, we can decouple and simplify the system considering *φ<sup>i</sup>* as an input, which provides a change exclusively in the *φ* output angle. Therefore, a nonlinear single-input single-output (SISO) system can be defined, having *φ<sup>i</sup>* inputs and *φ* outputs. Likewise, *θ<sup>i</sup>* and *δ<sup>i</sup>* inputs will affect only the output values of *θ* and *δ*, respectively, defining another two SISO systems.

Based on this, the soft neck can be modeled as three decoupled SISO systems. The transfer functions *Gθ*, *Gφ*, and *G<sup>δ</sup>* will model the actual outputs (*θ*,*φ*,*δ*) as a function of the new inputs (*θi*,*φi*,*δi*), defined by Equations (9)–(11). Given the simplifications we have considered, the real behavior will be different in several aspects, like showing interference between actuators and a nonlinear plant response. These effects will be discussed in the Experiments section.

Two different system identification methods are used. First, the set membership method as described in [32] is used for nonlinear identification, and second, recursive least-squares (RLS), as described in [33], is applied for different tilt configurations, which will result in a linear system for each RLS identification performed. The evolution of these systems, according to the inclination, will be studied.

In the case of RLS system identification, the new redefined inputs (*θi*, *φ<sup>i</sup>* and *li*) were considered instead of motor position inputs. Note that these are just the input redefinition, and the output angles still depend on the system dynamics. Although *f* functions are unknown, they are considered within the resulting models, although the nonlinear part may be neglected depending on the identification method.

## **3. Set Membership Non-Linear Identification**

This section briefly describes the Non-Linear Set Membership (NLSM) identification method proposed in [32].

Consider a system that has a Nonlinear AutoRegressive with eXogenous input (NARX) structure, as

$$y(k) = f\_o(\omega(k)) + \varepsilon(k) \tag{12}$$

where *ω*(*k*) is the system's regressor formed by past samples of the system inputs *u*1, *u*2 and the output *y*1, *y*2, as:

$$\omega(k) = [y\_i(k-1), \dots, y\_i(k-n\_{\underline{y}}), u\_i(k-1), \dots, u\_i(k-n\_u)]'\tag{13}$$

$$
\omega(k) \in \mathcal{W} \subseteq \mathcal{R}^n, n = \sum n\_{\mathcal{Y}\_i} + n\_{\mathcal{U}\_i} \tag{14}
$$

where *e*(*k*) represents the measurement noise and *W* is the function domain.

*i*

The NARX regressor is widely used in system identification considering its capacity of representing nonlinear dynamics and developing estimation algorithms which are computationally cost-efficient.

If *fo* is unknown, but a set of measurements of *yi*(*k*) and *ω*(*k*) are available for *k* = 1, ..., *N* and considering that the noise magnitude is bounded by :

$$|e(k)| \le \epsilon \tag{15}$$

and no statistical assumption on its behavior is made. The goal is to estimate ˆ *f* of *fo*, where ˆ *f* is the estimation of *f* .

Even though *f* is unknown, the following information is available:

$$f\_{\mathcal{O}} \in \mathcal{F} \doteq \left\{ f \in \mathbb{C}^1(\mathcal{W}) : \left\| f'(\omega) \right\| \leq \gamma, \forall \omega \in \mathcal{W} \right\} \tag{16}$$

where *f* (*ω*) denotes the gradient of *f*(*ω*) and *x* is the Euclidean norm. Therefore, we assume that the identified system is continuous on its first derivative and has maximum growth of *γ* for all the regressors applied to the function of interest.

On the other hand, if there is a Feasible System Set (FSS), which is the set of all systems in the space F which satisfies the following conditions:

$$FSS \doteq \left\{ \begin{array}{l} f \in \mathcal{F} : |y(k) - f(\omega(k))| \le \epsilon, \\ \text{and} \\ f \in \mathcal{F} : \frac{y(k) - y(k+1)}{\delta\_{\Gamma}} \le \gamma \end{array} \right\}, k = 1, 2, \ldots, N \tag{17}$$

therefore, there always exists a non-empty *FSS* and *fo* ∈ *FSS* when both assumptions on *fo* and *e* are true. Then, if we guarantee the validity of the conditions *γ* and over a set of measurements generated by the system to be identified, we will find a *FSS* = ∅. In [32], the procedure to guarantee conditions *γ* and over a data set is presented. For the following sections, prior assumptions are considered to be true.

Given that the aim of the model is to find the output generated by the system for a new input, it is necessary to distinguish the identification data set, *k*, and the new inputs *x*. Hence, for a given input *x* ∈ *W*, the optimal NLSM estimate of *fo*(*x*) is:

$$f\_c(\mathbf{x}) \doteq \frac{f\_{\rm ll}(\mathbf{x}) + f\_{\rm l}(\mathbf{x})}{2} \tag{18}$$

where:

$$f\_{\mathfrak{u}}(\mathbf{x}) = \min\_{1 \le k \le N} y(k) + \epsilon + \gamma \|\mathbf{x} - \omega(k)\|\tag{19}$$

$$f\_l(\mathbf{x}) = \max\_{1 \le k \le N} y(k) - \epsilon - \gamma \|\mathbf{x} - \omega(k)\|\tag{20}$$

As presented in [32],


$$f\_{opt} = \arg\inf\_{\hat{f}} \sup\_{f \in FSS} \left\| f - \hat{f} \right\|\_{p}.$$

The NLSM algorithm produces a non-linear, non-parametric model which is embedded on the data set. That is, there is no explicit equation that represents the input–output or physical variables relation.

For a new regressor value *<sup>x</sup>* ∈ R*n*, the NLSM model output *fc*(*x*) is evaluated through Algorithm 1.

In order to obtain the FSS, as described in [32,34], it is possible to execute the Algorithm 1 over the identification data set, updating the variable *γ* whenever the positive or negative projections *fu*(*i*), *fl*(*i*), over each data point *ω*(*i*) ∈ *ω*(*k*)∀*k* = *i* produces a greater, *fu*(*i*) < *y*(*i*), or lower, *fl*(*i*) > *y*(*i*), value.
