*3.1. Non-Linear Set Membership Data Set Generation*

In our specific problem, as described in Section 2, we have three motors that drive three tendons to actuate over the soft link and provide the desired pitch-roll motion. As described in Algorithm 1, to provide an estimation using Set Membership, we need to construct the FSS for a defined regressor that contains enough information so that all the identified system behaviors are contained in the FSS.

For this purpose, we define our regressor empirically since there are no standard methodologies to do so. Therefore, for simplicity, we run several system identifications using a neural networks MATLAB toolbox providing non-linear data-driven models for different regressor sizes. Once the NN is trained, we assume that the best neural network model uses the most informative regressor and corresponds to the original regressor selection. Later, to improve the computational time, the regressor is reduced by running different estimations modifying the number of elements in the regressor. In this way, less operations are required for each of the estimated data [35]. The chosen regressor is:

$$\begin{aligned} \omega(k) &= \quad [\underline{y}(k-1), \underline{u}\_1(k-2), \underline{u}\_1(k-3), \\ &\quad \underline{u}\_2(k-2), \underline{u}\_2(k-3), \\ &\quad \underline{u}3(k-2), \underline{u}3(k-3), \\ &\quad \underline{M}\_1(k-3), \\ &\quad \underline{M}\_2(k-3), \underline{M}\_3(k-3)] \end{aligned} \tag{21}$$

where *ui*(*k*) is the desired motor position for motor *i* at sample *k*, *Mi*(*K*) is the measured motor position at discrete time *k*, and *y* is the measured output at sample *k* when an estimator and control model is generated. Therefore, if noise or disturbances are detected in the measured output, the model aligns its dynamics using this information with the true model dynamics. On the other hand, if the required model is generated for prediction and simulation, the signal in (22) will be replaced by the previous model estimations at time *k* − *j y*(*k* − *j*). In this case, errors and disturbances detected at the output will not be perceived by the model unless a closed-loop control action modifies the input components of the regressor. In this case, if the model diverges from the real dynamics, they will not align with each other. The model architectures are given in Figure 3.

**Figure 3.** Model architectures. (**a**) An estimation and control model. (**b**) A prediction and simulation model.

With the regressor being defined as in (22), we generate our FSS by applying a sum of sinusoidals to each of the three motors such that the signals are not correlated and they give us a wide spectrum of the neck behavior. We capture the real motor positions, desired motor positions and, as output, the neck pitch and roll angles. Then, two separate models are generated.

The signal used for each motor *MP*(*i*) to create the FSS is described in Equation (22). The values used for the specific motor are listed in Table 2.

$$\begin{aligned} MP(i) &= \left( 0.6 \sin(2 \sin(\omega\_1 t + \phi\_1) + \cos(\omega\_1 t + \phi\_2)) \right) \* \dots \\ &\quad \* abs \left[ 3 + \sin(\omega 2t) + \sin(\omega\_3 t + \phi\_3) + \sin(\omega\_1 t - \phi\_4) + \dots \end{aligned} \tag{22}$$

$$\begin{aligned} \cos \dots &+ \sin(\omega\_4 t) + \sin(\omega\_5 t + \phi\_5) + \sin(\omega\_6 t + \phi\_6) \frac{\sin(\omega\_7 t + \phi\_7)}{2} \end{aligned} \tag{23}$$

As it can be seen, the different signals are non-linear and have different frequencies and phases. This provides a pseudo-random interaction and covers a wide range of operational modes of the neck. The frequency spectrum that was chosen is coherent with the system bandwidth, which is 4 (rad/s), generating soft, continuous, and human-like motion in the axes of interest. By modifying the phase and mixing the minimum 0.25 (rad/s) and maximum rad/s frequencies, we aim to explore in a single experiment a wide range of motions providing sufficient information for the FSS. However, it is necessary to point out that some system dynamics did not occur during the proposed study scenario, such as the three tendons pulling at the same time with the same force, which keeps a static neck position with different stiffness, as well as continuous single tendon activation, to name some. Even if the proposed FSS does not cover the full system dynamics, the chosen signal should cover the normal operational range for the soft neck.


**Table 2.** Values used for the identification data set creation.

The identification and validation data sets are given in Figure 4, where 10,000 samples were taken, 7000 for the FSS (Identification Data Set) and 3000 for the validation set.

As seen in the Figure, the rotational position of the tendon is maximum 6 rad. Having a pulley diameter of 15 mm, each tendon has linear displacement of the tendon 7.5 mm/rad, and therefore, a maximum linear displacement of ≈45 mm. This generates inclinations of ±20 deg for the pitch and [−20, 40] deg for the roll. This provides a wide dynamic range of motion. In addition, the motion frequency was set to replicate a human-like motion which is the region of interest.

**Figure 4.** Input and output signals.

#### **4. Recursive Least-Squares Linear Identification**

Despite the nonlinear nature of the plant, a simple recursive least-squares (RLS) identification was performed. On the one hand, it will show a qualitative estimation of the system's nonlinearity degree, and on the other hand, a linear model could be used to solve the nonlinearity issues by means of robust or adaptive control strategies.

The RLS identification algorithm can be described using an ARX structure:

$$\hat{y}(t) = -a\_1 y(t-1) - \dots - a\_{n\lambda} y(t-na) + b\_1 u(t-1) + \dots + b\_{nb} u(t-nb),\tag{23}$$

with *y*(*t*) and *u*(*t*) being the plant output and input variables, with a matrix representation as follows:

$$
\hat{y}(t) = \theta \phi'(t-1) \tag{24}
$$

$$\theta = [a\_{1\prime} \dots a\_{na\prime} b\_{1\prime} \dots b\_{nb}] \tag{25}$$

$$\phi(t-1) = \left[ -y(t-1), \dots, -y(t-na), u(t-1), \dots, u(t-nb) \right] \tag{26}$$

Increasing one time-index (*y*ˆ(*t* + 1) = *θφ* (*t*)), Equations (24)–(26) provide the output prediction, based on the model parameters (*θ*), and the set of past inputs and outputs (*φ*(*t* − 1)). Comparing the next actual system output with this predicted value results in the prediction error:

$$
\epsilon(t) = y(t) - \hat{y}(t) \Rightarrow \epsilon(t+1) = y(t+1) - \hat{y}(t+1) \tag{27}
$$

In order to minimize this error, different algorithms can be used. In the least-squares case, the squared sum of all errors is the variable to be minimized. Since the parameters that minimize the error produced by the least-squares solution can also be obtained from the preceding parameters (recursively), the algorithm can be expressed using recursion, as shown below:

$$
\hat{\theta}(t+1) = \hat{\theta}(t) + F(t+1)\phi(t)\epsilon(t+1) \tag{28}
$$

$$F(t+1) = F(t) - \frac{F(t)\phi'(t)\phi(t)F(t)}{1 + \phi(t)F(t)\phi'(t)}\tag{29}$$

$$
\epsilon(t+1) = y(t+1) - \theta(t)\phi'(t) \tag{30}
$$

These equations define the operations needed to find ˆ *θ*(*t* + 1) based on the previous parameters and captured inputs and outputs. An improved method is described in [33] as RLS with a constant forgetting factor (CFF-RLS). It will not be necessary at this point, as the soft neck system identification will be done offline, but it can be used in the future, if an adaptive scheme is proposed as a solution to the nonlinearity issues. See [33] or [36] for a more detailed discussion about RLS and other identification methods.

Using the described inputs and outputs definition of Section 2, the same data captured in the identification experiments were used in order to obtain a plant model. As the data capture is based on the motor positions, we can find the equivalent input values using Equations (9) and (10), and consider these inputs. The outputs will be the same as in the other cases, the neck pitch and roll angles.

For example, Figure 5 shows part of the inputs and outputs considered for the pitch and roll RLS identification.

**Figure 5.** Examples of decoupled identification datasets. Redefined input targets compared to measured angles for pitch (**left**), and roll (**right**).

A small delay can be observed in the data set example shown in that Figure. This delay must be considered during the system identification process. The resulting transfer functions corresponding to the pitch dynamics, *G<sup>θ</sup>* (*s*), and roll dynamics, *Gφ*(*s*), obtained using the RLS algorithm through the entire data set, are:

$$\mathbf{G}\_{\theta}(\mathbf{s}) = \mathbf{e}^{-0.08\epsilon} \frac{27.691}{(\mathbf{s} + 5.293)} \quad \mathbf{G}\_{\phi}(\mathbf{s}) = \mathbf{e}^{-0.08\epsilon} \frac{25.424}{(\mathbf{s} + 4.938)}.\tag{31}$$

The model unit input time responses are shown in Figure 6 for the described system model. Note how both systems' (*G<sup>θ</sup>* (*s*), *Gφ*(*s*)) static gains are close to 5, showing a stationary response above the unit input level, as expected from Figure 5.

**Figure 6.** Unit input time response for *G<sup>θ</sup>* (*s*) (**left**) and *Gφ*(*s*) (**right**).

Once the linear models of the soft neck decoupled system are determined, a new input response simulation was performed using those models, together with a new data set for validation and accuracy check. A partial plot of these results is shown in Figure 7 for the pitch and roll angles.

**Figure 7.** Validation example of the identified models for *G<sup>θ</sup>* (*s*) (**left**) and *Gφ*(*s*) (**right**).

Note how although the linear model captures the system behavior quite well, there are mismatches due to plant non-linearity. In order to deal with these problems, a robust controller could be used in the future, since it will provide a constant behavior despite plant parameter changes or non-linearities.

### **5. Experimental Setup**

The main objective of this work was to extract a 3D model for the soft neck platform. The selected identifications were done in an offline configuration, with an open-loop data capture scheme. In this sense, a set of experiments has been designed for data capture.

As stated before, as we expected a correlation between the motor position and the inclination of the top platform, the motors' states are therefore considered as inputs, while the measured sensor angles are used as system outputs. The captured data involve the following inputs and outputs:


Neck actuator motion was programmed to follow a composition of sinusoidal functions, as described in Section 3.1. The captured motion describes human-like movement. Input and output sets can be seen in Figure 4.

All models used the same data for modeling the system, enabling direct comparisons of the models. Additional tests were also captured for validation purposes.
