Constrained methods are widely diffused, for example, in the case of offline trajectory planning. The general approach consists in the formulation of a constrained optimization problem that should optimize a given objective, such as minimum execution time [1] or minimum energy consumption [2]. This is also the case of global redundancy resolution methods [3], where the robot has to perform an assigned task, and the redundancy can be exploited to optimize a desired objective, such as maximum manipulability [4], dexterity [5], or maximum joint range availability [6]. However, all these methods are typically performed offline as they need the prior knowledge of the task to be performed, and they lead to heavy computational burdens that do not permit their online execution. For this reason, they are not able to handle online trajectory generation and re-planning techniques and this is a relevant limitation, considering recent robotic applications such as collaborative robotics and robots in unstructured environments. To handle the robot constraints in online planning and control methods,

optimization-based techniques are becoming more and more widespread as the increasing computing power of modern processors allows their real-time implementation with small sampling periods.

Quadratic programming, for example, is intensively exploited in the resolution of the Inverse Kinematic (IK) problem of redundant manipulators [7–11]. In brief, a manipulator is termed kinematically redundant with respect to a task when it has a number of degrees of freedom greater than the one of the desired task. The IK problem does not present a unique solution, as the transformation between the joint and the task space results in an undetermined system. Such problem is usually addressed at differential level in order to exploit the linear relation between joint and task velocities given by the Jacobian of the robot. This can be formulated as a chain of Quadratic Programs (QP) with decreasing priorities: the tasks are written as cost functions and the robot limits are represented as linear constraints. The possibility of including robot constraints in the QP represents a great advantage with respect to other typical methods (such as the unconstrained projection of secondary tasks in the null space of the Jacobian [12,13]). Further developments in this field explicitly face the problem of possible saturations and task deformations by introducing a scaling variable into the QP to slow down the execution of the tasks if needed [14].

A further improvement of these methods is represented by the use of Model Predictive Control (MPC) techniques, which aim at overcoming issues related to the fact that the previous methods only consider the current state of the robot. In fact, MPC-based methods are able to take into account the future evolution of the system, the tasks, and the constraints, improving the behavior of the robot in terms of task satisfaction and smoothness [15,16]. MPC is also applied to constrained motion planning [17,18] or low-level control applications [19].

Disregarding the particular field of application, a fundamental issue to address when using optimization-based control techniques is the feasibility of the online optimization problem. In fact, feasibility must be ensured in any possible state of the system. Otherwise, a solution to the problem might not exist, and this could cause the algorithm to stop (unless specific but sub-optimal strategies are activated when infeasibility occurs).

Feasibility of the online optimization problem is strictly related to the concept of set viability [20]. A set is said to be viable if, given an initial state within such set, the state trajectory can be kept within the set by means of a proper and realizable input function. In other words, keeping the state within a viable set ensures feasibility for all future time instants. On the contrary, if the state exits the viable set, infeasibility of the control problem will surely occur at a certain time in the future. This is particularly relevant in the case of state-constrained systems, as infeasibility might occur if no viability conditions are added to the control problem but only state and input limits are taken into account.

For robot manipulator control, this issue is common to any control strategy that takes into account both position and acceleration/torque constraints. In this case, the simplistic imposition of box constraints on the position and the acceleration/torque of the robot does not ensure the existence of a feasible solution to the optimization problem. For example, if a robot joint approaches its position limit with high velocity, the position bound will be exceeded due to the bounded admissible deceleration. Few works addressed this issue by means of manually-tuned heuristic strategies that aimed at reducing the velocity of the robot when it approached its position limits [21–23]. Many other control strategies did not take into account position bounds, assuming that the reference trajectories were implicitly feasible [24,25]. Notice that the viability guarantee does not ensure the feasibility of the reference trajectory, which could be only verified offline when the whole trajectory is known a priori. Such approach does not apply to online methods, which should be able to handle online trajectory generation and re-planning. In these cases, we can only ensure the feasibility of the optimal control control, which typically aims at minimizing the deviation of the performed motion with respect to the reference trajectory. This means that, although the desired motion might be infeasible, the optimal control problem remains feasible. In this case, a deformation with respect to the nominal trajectory could arise.

A formal viability guarantee was given in [26], which proposed an invariance control scheme for constrained robot manipulators. In particular, viability conditions were derived in order not to exceed workspace limitations with bounded joint accelerations. In [14], a local redundancy resolution technique was proposed: to respect a given deceleration limit in the breaking phase, the position-velocity state region was limited by imposing that the robot was able to stop within the given joint limit (in fact, limiting the state region to a viable set). This approach was similarly adopted in [27] in the context of robust constrained motion planning. Recently, viability for joint-constrained manipulators was explicitly addressed in [28], with a particular focus on the discrete-time implementation of the viability constraints.

All the above-mentioned approaches were based on the derivation of an analytical viability condition for a double integral system with bounded input and output. The resulting condition is quadratic in the velocity state. However, as all the above-mentioned works refer to local/feedback methods, such condition could be easily linearized around the current system state, resulting in variable box constraints on the control actions. However, such linearization is not possible for predictive strategies. In this case, a linear approximation of the viability set should be obtained. For example, Faroni et al.[29] gives a sufficient condition to approximate the quadratic viable set with a linear constraint and maintaining the viability property. Although the approximation technique is easy to implement and does not remarkably increase the computational complexity of the problem, the resulting allowed state region could result to be significantly smaller than the original one.

In this paper, we propose a method to derive a viable convex polyhedron for a robotic system with bounded joint positions, velocities, and accelerations. In particular, a simple optimization problem is set up to determine the maximum polyhedron that approximates the original viable set. The original quadratic condition is approximated by a polyhedron with a given number of sides, by maximizing the area of the allowed velocity-position state-region. Moreover, viability of the resulting set is ensured by imposing that, for all points on the polyhedron boundary, there exists a realizable input action that keeps the next state within the polyhedron itself, which is also demonstrated to be convex and can be therefore rewritten as a linear inequality constraint in optimization-based algorithms (such as linear MPC techniques). The paper shows that, by increasing the number of sides, the polyhedron gives a better approximation of the original viable set, as expected. This means that a larger admissible state region can be exploited by the controller when the resulting constraints are included in the optimization problem. Consequently, the controller can obtain smaller tracking errors when the desired task requires the robot joint to get close to the maximal viable set boundary. In particular, numerical results show the enlargement of the admissible state region as the number of sides increases. Such improvement is more and more evident as the maximum acceleration values gets smaller. Finally, experimental results on a Universal Robots UR10 manipulator demonstrates the validity of the proposed approach. In particular, an MPC algorithm is applied to the tracking problem of a given joint-space trajectory and the results with different viability conditions are compared. Firstly, an experimental example shows that the use of an invariance condition is of vital importance to ensure both the feasibility of the online optimization problem and the satisfaction of the manipulator limits. Secondly, the results show that the viable sets obtained by means of the proposed method permits to obtain smaller (or null) tracking errors in the case the required motion is close to the viability boundary, increasing the performance of the MPC controller.

The paper is organized as follows. Section 2 introduces the concept of set viability applied to constrained kinematic control of robot manipulators. Section 3 illustrates the proposed method for the computation of the optimal viable polyhedron and gives viability and convexity proves. Numerical results demonstrate the effectiveness of the proposed approach in Section 4, while experimental results on a Universal Robot UR10 manipulator are shown in Section 5. Section 6 concludes the paper.

#### **2. Feasibility of the Constrained Kinematic Control Problem**

Consider a generic robot joint and denote with *q*, *q*˙, *q*¨ its position, velocity, and acceleration, respectively. The limits on *q*, *q*˙, and *q*¨ can be therefore expressed as:

$$\mathfrak{q}\_{\text{min}} \stackrel{\kappa^\*}{=} \mathfrak{q} \quad \stackrel{\kappa^\*}{\leq} \mathfrak{q}\_{\text{max}}.\tag{1}$$

$$
\dot{q}\_{\text{min}} \le \dot{q} \le \dot{q}\_{\text{max}}.\tag{2}
$$

$$
\vec{q}\_{\text{min}} \le \vec{q} \le \vec{q}\_{\text{max}}.\tag{3}
$$

where *q*min, *q*˙min, *q*¨min, and *q*max, *q*˙max, *q*¨max are the minimum and maximum joint position, velocity, and acceleration, respectively.

Assume a discrete-time implementation with sampling period *T* in which the acceleration is considered as constant along each sampling period. At time *k*, kinematic limits at the next sampling time *k* + 1 can be easily written as linear inequalities in the joint acceleration *q*¨*k* as:

$$q\_{\rm min} \le \left| q\_k + T \, \dot{q}\_k + \frac{1}{2} \, T^2 \, \ddot{q}\_k \right| \le q\_{\rm max} \tag{4}$$

$$
\dot{q}\_{\min} \le \left. \dot{q}\_k + T \, \dot{q}\_k \right| \le \dot{q}\_{\max} \tag{5}
$$

$$
\not\models \not\models\_{\text{min}} \leq^\* \not\models \not\models\_{\text{max}} \cdot \text{ } \tag{6}
$$

However, the simplistic imposition of Equations (4)–(6) may result in an empty admissible set, as no feasible solution might exist that satisfies all constraints at the same time [28].

Indeed, the existence of a solution is guaranteed if Equations (4)–(6) are feasible for all the admissible states of the system. This is strictly correlated to the concept of set viability, as mentioned in Section 1. As each joint is modeled as a double integrator, the viability analysis traces back to the viability of a double integrator system with bounded input and output. By imposing the feasibility of Equations (4)–(6), the maximal viability set for the double integrator can be derived analytically [26]. Intuitively, it can be calculated by imposing that, applying the maximum deceleration *q*¨min, the joint stops at *q* = *q*max with null velocity (*q*˙ = 0). The resulting condition (valid for the upper bound) is given by:

$$\begin{cases} q - \frac{q^2}{2q\_{\text{min}}} - q\_{\text{max}} \le 0 & \text{if } \dot{q} > 0\\ q - q\_{\text{max}} \le 0 & \text{otherwise} \end{cases} \tag{7}$$

which expresses a quadratic condition in the system states. An analogous condition for the lower bound can be derived likewise.

Such condition can be easily linearized around the current velocity and position, in the case of local methods, as in [14,27]. The inclusion of Equation (7) in the QP ensures that the states of the system (i.e., the joint velocities and positions) remain within a viable set for which the problem is feasible. More details about the discrete implementation in robotic systems can be found in [28].

The necessity of deriving a linear approximation of the original quadratic equation (Equation (7)) comes from the advantages given by the linear formulation in optimization-based controllers (in particular, the significant decrement of computational time and the ease of implementation of linear MPC with respect to the nonlinear one). However, in the case of predictive strategies, the linearization adopted for local methods is not applicable.

To tackle this issue, Faroni et al. [29] proposed a linear viability condition based on a single constraint for each joint. Firstly, it showed that approximating the quadratic constraint in Equation (7) by means of a straight line passing through the extreme points of the maximal viable set (i.e., (*q*max, 0) and the intersection between Equation (7) and *q*˙ = *q*˙max in Figure 1) leads to a non-viable set. Then, it derives a linear viability condition by imposing that the maximum deceleration along the linear

constraint is exactly equal to the maximum admissible deceleration. In this way, the joint states can be always kept within the resulting set by means of a realizable acceleration. Such condition is given by:

$$-\frac{\ddot{q}\_{\text{min}}}{\dot{q}\_{\text{max}}}q+\dot{q}+\frac{\ddot{q}\_{\text{min}}}{\dot{q}\_{\text{max}}}q\_{\text{max}}\leq 0\tag{8}$$

for the upper bound, and likewise for the lower one.

Inequalities (Equations (7) and (8)) are graphically represented in Figure 1. Linearity of Equation (8) is a great advantage, as it allows the straightforward inclusion in an optimization problem as a linear constraints. However, it is clear that Equation (8) implies a conservative reduction of the available state-space region compared to Equation (7), and such reduction worsens as the acceleration limits become smaller.

From a practical perspective, the reduction of the admissible state region shows its drawbacks when the robot is required to perform a motion that would violate Equation (8). In fact, the states laying between the linear and the quadratic constraints in Equations (7) and (8) are potentially realizable by the robot, but they are automatically excluded by the controller, with consequent deformation of the desired trajectory. Similarly, for redundant manipulators, shrinking the admissible state region could results in a worse satisfaction of the secondary objectives.

**Figure 1.** Viable admissible set for the double integrator system. Black: Maximal region, given by the quadratic inequality in Equation (7). Blue: Conservative linear inequality in Equation (8), as proposed in [29].

#### **3. Proposed Method**

As mentioned in the previous section, viability conditions are fundamental to ensure the feasibility of the control problem. The derivation of linear viability conditions (e.g., Equation (8)) allows their straightforward implementation in an optimization-based control framework, such as MPC algorithms. A strategy to enlarge the admissible region without waving its linearity is therefore proposed. It consists in increasing the number of linear constraints that compose the admissible set. We can therefore state the following problem.

**Problem 1.** *Given the acceleration, velocity, and configuration limits for each joint, find the maximal viable convex polyhedron determined by a given number of linear inequalities.*

We address this problem by converting it into an optimization problem. In particular, for each joint, two optimization problems are set up (one for the upper position bound and one for the lower position bound). Without loss of generality, consider the upper configuration bound for a generic joint (i.e., the first quarter in Figure 1). The optimization variables of the problem are given by the coordinates of the extremes of the segment composing the polyhedron shown in Figure 2. Given a number of segments *h*, such variables are denoted with the vectors:

$$p\_{\mathbf{x}} = (p\_{\mathbf{x},0}, \dots, p\_{\mathbf{x},\mathbf{j}}, \dots, p\_{\mathbf{x},\mathbf{h}}),\tag{9}$$

$$p\_y = (p\_{y,0}, \dots, p\_{y,j}, \dots, p\_{y,h}),\tag{10}$$

whereas a generic point on the polyhedron is denoted with *Pj*(*px*,*j*, *py*,*j*).

Now, impose that (*px*,0, *py*,0)=(*q*mean, *q*˙max), where *q*mean = *q*max+*q*min 2 and (*px*,*h*, *py*,*h*)=(*q*max, 0). To maximize the area covered by the polyhedron, a cost function equivalent to the opposite of such area is defined as follows:

$$\psi = \frac{1}{2} \sum\_{j=1}^{l} (p\_{x,j-1} - p\_{x,j})(p\_{y,j-1} + p\_{y,j}).\tag{11}$$

**Figure 2.** Construction of the polyhedron P1.

The extremes of the segments are required to lay in the first quarter and have to respect the maximum position and velocity bounds. This results in the following box constraints:

$$\mathfrak{gl} \text{ means } \stackrel{\kappa}{\le} \mathfrak{p} x \stackrel{\kappa}{\le} \mathfrak{gl} \text{max}\_{\mathsf{M}} \tag{12}$$

$$0 \le p\_{\underline{y}} \le \dot{q}\_{\text{max}}.\tag{13}$$

Moreover, as a first necessary condition to the convexity of the polyhedron, we impose:

$$p\_{\mathbf{x},j} - p\_{\mathbf{x},j-1} \ge 0 \quad \forall j = 1, \dots, h,\tag{14}$$

$$p\_{y,j-1} - p\_{y,j} \ge 0 \quad \forall j = 1, \ldots, h. \tag{15}$$

One last constraint needs to be imposed to ensure that the resulting polyhedron is viable and convex. To this purpose, we resort to a geometrical reasoning. Consider a generic segment *Pj Pj*+1 whose extremes are the points (*px*,*j*, *py*,*j*) and (*px*,*j*+1, *py*,*j*+1) and denote with *nj* the unitary vector normal to the segment and directed toward the inner region of the polyhedron. It results that:

$$\vec{m}\_{\vec{j}} = \left( \frac{p\_{y,\vec{j}+1} - p\_{y,\vec{j}}}{||\overline{P\_{\vec{j}}P\_{\vec{j}+1}}||}, \frac{p\_{x,\vec{j}} - p\_{x,\vec{j}+1}}{||\overline{P\_{\vec{j}}P\_{\vec{j}+1}}||} \right). \tag{16}$$

*Robotics* **2018**, *7*, 41

Consider then a generic point *Pζ* (*px*,*ζ* , *py*,*ζ* ) on the segment. By applying the maximum deceleration *q*¨min to the system, the states are forced to move on the curve given by:

$$\begin{cases} q = p\_{x,\tilde{\Box}} + p\_{y,\tilde{\Box}}t + \frac{1}{2}\ddot{q}\_{\text{min}}t^2\\ \dot{q} = p\_{y,\tilde{\Box}} + \ddot{q}\_{\text{min}}t \end{cases} \tag{17}$$

The tangent vector to such curve in the (*q*, *q*˙)-plane is given by:

$$\vec{\sigma} = \left(\frac{dq}{dt}, \frac{d\vec{q}}{dt}\right) = \left(p\_{y,\vec{\zeta}} + \vec{q}\_{\text{min}}t, \,\vec{q}\_{\text{min}}\right). \tag{18}$$

The vector*t* at *t* = 0 is therefore the tangent vector in the point *Pζ* and results to be:

$$
\vec{t}(p\_{x,\mathbb{Z}\_{\mathcal{I}}}, p\_{y,\mathbb{Z}\_{\mathcal{I}}}) = (p\_{y,\mathbb{Z}\_{\mathcal{I}}}, \vec{q}\_{\min}).\tag{19}
$$

Note that*t*(*px*,*ζ* , *py*,*ζ* ) represents the direction of the state movement when the maximum realizable deceleration is applied. The state will therefore be able to remain below the considered segment if the following condition holds:

$$
\vec{n}\_{\rangle} \cdot \vec{t}(p\_{x,\zeta}, p\_{y,\zeta}) \ge 0 \qquad \forall \, (p\_{x,\zeta}, p\_{y,\zeta}) \in \overline{P\_{\rangle} \overline{P\_{\rangle+1}}} \tag{20}
$$

where (·) denotes the scalar product between two vectors. As

$$\min\_{p\_{x,\zeta}, p\_{y,\zeta}} \left( \vec{n}\_{\rangle} \cdot \vec{t}(p\_{x,\zeta}, p\_{y,\zeta}) \right) = \vec{n}\_{\rangle} \cdot \vec{t}(p\_{x,\zeta}, p\_{y,\zeta}),\tag{21}$$

condition in Equation (20) can be imposed on the sole extreme point *Pj* of each segment and becomes:

$$p\_{y,j}(p\_{y,j+1} - p\_{y,j}) + \ddot{q}\_{\text{min}}(p\_{x,j} - p\_{x,j+1}) \ge 0. \tag{22}$$

Finally, the resulting optimization problem can be written as:

$$\begin{aligned} \underset{p\_x, p\_y}{\text{minimize}} \quad & \forall \\ \text{subject to} \quad & q\_{\text{mean}} \le p\_x \le q\_{\text{max}} \\ & 0 \le p\_y \le q\_{\text{max}} \\ & p\_{x,j} - p\_{x,j-1} \ge 0 \quad \forall j = 1, \dots, h \\ & p\_{y,j-1} - p\_{y,j} \ge 0 \quad \forall j = 1, \dots, h \\ & p\_{y,j}(p\_{y,j+1} - p\_{y,j}) + \bar{q}\_{\text{min}}(p\_{x,j} - p\_{x,j+1}) \ge 0 \quad \forall j = 1, \dots, h - 1 \\ & p\_{x,0} = q\_{\text{mean}} \\ & p\_{y,0} = q\_{\text{max}} \\ & p\_{x,h} = q\_{\text{max}} \\ & p\_{y,h} = 0 \end{aligned} \tag{23}$$

*Robotics* **2018**, *7*, 41

Following the same reasoning, an analogous optimization problem can be set up for the lower position bound. Note that, if the joint limits are symmetric, the solution for the lower bounds can be obtained by "mirroring" the solution of Equation (23) into the third quarter of the (*q*, *q*˙)-plane. Notice that, as Equation (23) is non-convex, global optimization such as genetic or multi-start algorithms should be adopted to solve the optimization problem (see Section 4).

Denoting with P1 the polyhedron determined by the segments obtained as solution of Equation (23), and with P2 the analogous polyhedron obtained for the lower bounds, the overall polyhedron for the single joint is given by:

$$\mathcal{P} = \mathcal{P}\_1 \cup \mathcal{P}\_2 \cup \ [q\_{\text{min}}, q\_{\text{mean}}] \times [0, \dot{q}\_{\text{max}}] \cup \ [q\_{\text{mean}}, q\_{\text{max}}] \times [\dot{q}\_{\text{min}}, 0]. \tag{24}$$

**Proposition 1.** *Polyhedron* P *is a viable set with respect to the given joint position, velocity and acceleration limits. Moreover, such polyhedron is convex.*

**Proof (Viability).** Viability of the set is implicitly ensured by the imposition of Equation (22).

**Proof (Convexity).** To prove it by contradiction, denote the solution of Equation (23) with *p*∗ *x*, *p*∗ *y* and assume that the corresponding polyhedron P∗ 1 (i.e., the part of the polyhedron in the first quarter) is non-convex. Consider then a triple of consequent extremes {*Pj*−1, *Pj*, *Pj*+1}, which causes a non-convexity. Equations (20) and (21) imply:

$$\left(\vec{n}\_{j-1} \cdot \vec{t}(p\_{x,j-1}^\*, p\_{y,j-1}^\*) \ge 0,\tag{25}$$

$$\left(\vec{x}\_{j} \cdot \vec{t}(p\_{x,j}^{\*}, p\_{y,j}^{\*})\right) \geq 0. \tag{26}$$

Consider now the point *Pc* given by the projection of *Pj* onto *Pj*−1 *Pj*+1 and denote with *nc* the vector normal to *Pj*−1 *Pj*+1. We want to prove that the polyhedron obtained by substituting *Pc* to *Pj* in P∗ 1 is a feasible and more efficient solution of Equation (23). First, constraints in Equations (12)–(15) are straightforwardly satisfied. Moreover,

$$\vec{n}\_{\varepsilon} \cdot \vec{t}(p\_{x,j-1}^\*, p\_{y,j-1}^\*) \ge \vec{n}\_{j-1} \cdot \vec{t}(p\_{x,j-1}^\*, p\_{y,j-1}^\*) \ge 0 \tag{27}$$

as *nj*−1 ·*i* ≤ *nc* ·*i* ≤ 0 and*nc* ·*j* ≤ *nj*−1 ·*j* ≤ 0 by construction (where*i* and*j* denote the unitary vectors directed as the horizontal and the vertical axis, respectively).

Moreover, as *q*¨min in Equation (19) is assumed to be negative, the following inequality holds:

$$\|\vec{m}\_{\varepsilon} \cdot \vec{\mathbf{f}}(p\_{x,\varepsilon}, p\_{y,\varepsilon})\| \ge \vec{m}\_{\varepsilon} \cdot \vec{\mathbf{f}}(p\_{x,j-1}^\*, p\_{y,j-1}^\*) \ge 0. \tag{28}$$

Note that Equations (27) and (28) implies that the new candidate solution satisfies Equation (22) and is therefore a feasible solution to Equation (23). Moreover, the candidate solution is more efficient than the assumed one, as

$$\Psi(P\_{0\prime},\ldots,P\_{j-1},P\_{\mathbf{c}\prime},P\_{j+1},\ldots,P\_{\mathbf{h}}) \le \Psi(P\_{0\prime},\ldots,P\_{j-1},P\_{j\prime},P\_{j+1},\ldots,P\_{\mathbf{h}}),\tag{29}$$

which means that *p*∗ *x*, *p*∗ *y* are not the optimal solution of Equation (23), as a feasible and more efficient solution that eliminates such non-convexity exists. Applying this reasoning to any triples of extremes that cause a non-convexity implies that the optimal polyhedron P1 must be convex. The same demonstration can be applied to P2 leading to the same conclusion.

#### **4. Numerical Results**

Consider the first joint of the robot manipulator Universal Robot UR10. The joint configuration and velocity limits from the datasheet are given by:

$$q\_{\text{max}} = -q\_{\text{min}} = 3.14 \text{ rad}, \qquad \dot{q}\_{\text{max}} = -\dot{q}\_{\text{min}} = 2.16 \text{ rad/s}. \tag{30}$$

Typically, acceleration limits are not explicitly given in the datasheet. Indeed, the actual dynamic limit of the joint is usually determined by the available joint torque. However, handling viability directly in terms of torque limits is typically avoided due to the complexity of the problem. The common approach consists in estimating a (usually conservative) corresponding acceleration limit for each joint separately. A method for the estimation of such acceleration limit is given, for instance, in [28]. As a first example, assume that the acceleration limits are given by *q*¨max = −*q*¨min = 6 rad/s2.

We solve the optimization problem (Equation (23)) for different values of *h* to evaluate the resulting polyhedrons as the number of edges grows and to compare the results to the quadratic and the linear viability constraints (Equations (7) and (8), respectively). Table 1 shows the values of the inner area of polyhedron P1. The values are normalized with respect to the maximal viable set obtained by means of Equation (7), whose area is therefore equal to one (last column of Table 1). As expected, the smaller area is given by Equation (8) (first column of Table 1), while the extension of the polyhedron obtained by solving Equation (23) grows with the number of edges. In other words, the larger the number of edges, the closer the optimal polyhedron is to the maximal extension (given by Equation (7)). This is clearly shown in Figure 3, which depicts the different viable sets.

An analogous example is performed by using *q*¨max = −*q*¨min = 3 rad/s2. The values of the normalized area of P1 are shown in Table 2 and the resulting viable sets are depicted in Figure 4. The results lead to conclusions similar to the ones given by the previous example. However, it is clear that the magnitude of the phenomenon grows as the acceleration limits get smaller.

Of course the use of a larger value of *h* also gives some drawbacks. First, the computational complexity of Equation (23) rapidly grows with *h* as the number of variables and constraints is linear in *h*. As an example, the time needed to solve Equation (23) for *h* = 1 and *h* = 11 was in the order of 1 second and 15 s, respectively (the computation was performed in Matlab using a multi-start gradient-descent method on a standard laptop mounting a 2.5 GHz Intel Core i5-2520M processor).

Furthermore, the value of *h* determines the number of linear inequalities describing the polyhedron. In optimization-based control algorithms (such as linear MPC), this affects the computational complexity of the online optimal control problem, which is typically a critical issue in online methods, as mentioned in Section 1.

**Table 1.** Inner area of P1 for different values of *h* in case *q*¨max = −*q*¨min = 6 rad/s2.

0.9139 0.9636 0.9773 0.9836 0.9871 0.9895 0.9911 0.9920 0.9923 0.9925 1

**Figure 3.** Comparison of the viable polyhedrons for different values of *h* (portion of interest of the first quarter) in the case *q*¨max = −*q*¨min = 6 rad/s2.

**Figure 4.** Comparison of the viable polyhedrons for different values of *h* (portion of interest of the first quarter) in the case *q*¨max = −*q*¨min = 3 rad/s2.

#### **5. Experimental Results**

Experimental were also performed on a six-degrees-of-freedom Universal Robot UR10 manipulator to prove the validity of the proposed approach. For the purposes of the paper, the joint position and velocity limits have been set to slightly conservative values with respect to the ones given in the datasheet.  The experimental platform is shown in Figure 5. The objective of these experiments was the evaluation of the performance of an optimization-based control algorithm when different viability constraints are implemented. To this purpose, a simple MPC scheme was applied to a trajectory following problem. The MPC control scheme is in charge of following a position reference signal in the joint space. The model implemented in the MPC consists of a double integrator for each joint, as typical of robot kinematic control [15,30]. The input action is therefore represented by the joint accelerations, which then feed the low-level controller of the robot. Joint position, velocity, and acceleration limits can be implemented as linear constraints in the MPC (see [15] for details) and the online optimal control problem results to be a QP which minimizes the weighted sum of the tracking error and the control effort, as typical of linear MPC [31]. Notice that the choice of such a simple control scheme is due to the will of highlighting the behavior of a linear MPC controller in the case of different viability conditions. However, the proposed method could be straightforwardly applied to more complex MPC techniques such as [15,16,29].

**Figure 5.** Universal Robot UR10 used as experimental platform.

The MPC algorithm was implemented using a sampling period *T* = 8 ms and by setting a predictive and a control horizon *N* = 20 sampling periods. The robot trajectory is controlled by means of a ROS-based control architecture. Namely, a position controller runs in a ROS Kinetic Ubuntu 16.04. The controller communicates with the robot by means of a TCP connection. The controller takes the MPC position output as reference and receives the actual joint position. The controller output is the sum of a proportional action, with gain equal to seven, and a feedforward term equal to the MPC velocity output.

*Robotics* **2018**, *7*, 41

We considered a simple trajectory given by a straight line in the joint space, parameterized with respect to the normalized longitudinal length along the path, denoted with *r* ∈ [0, 1]. The trajectory *Q*des(*t*) is therefore defined as:

$$Q\_{\rm des}(t) = Q\_{\rm start} + \left(Q\_{\rm end} - Q\_{\rm start}\right)r(t),\tag{33}$$

where *Q*start = ( 0, −1, 0.5, 0, 0, 0.5 ) rad and *Q*end = ( 2.9, −0.5, −1.5, 1, −0.5, 2.9 ) rad are the initial and the final points of the trajectory, respectively, and *r* : [0, *t*end] → [0, 1], *t* → *r*(*t*) is defined as a timing law with trapezoidal velocity profile, where *t*end represents the total time of the trajectory, as depicted in Figure 6.

**Figure 6.** Timing law with trapezoidal velocity profile.

The control problem the MPC controller has to solve at each cycle *k* therefore results:

$$\begin{aligned} \underset{Q}{\text{minimize}} & \quad \sum\_{i=1}^{N} \left\| Q\_{\text{des}}(k+j) - Q(k+j) \right\|^2 + \lambda \sum\_{i=0}^{N-1} \left\| \vec{Q}(k+j) \right\|^2\\ \text{subject to} & \quad Q\_{\text{min}} \le Q(k+j) \le Q\_{\text{max}} \qquad \forall j = 1, \ldots, N\\ & \quad Q\_{\text{min}} \le Q(k+j) \le Q\_{\text{max}} \qquad \forall j = 0, \ldots, N-1\\ & \quad Q(k+1) = Q(k) + T\dot{Q}(k) + 0.5T^2 \dot{Q}(k)\\ & \quad \dot{Q}(k+1) = \dot{Q}(k) + T\ddot{Q}(k) \end{aligned} \tag{34}$$

where *Q*, *Q*˙ , and *Q*¨ ∈ R6 are the joint position, velocity, and acceleration vectors, respectively, and the controller effort weighting factor *λ* = 1 × 10−6 was tuned empirically. Different implementations of the position limits will then be added to the problem in order to evaluate the different behaviors of the system.

As a first example, the total time of the trajectory was chosen as *t*end = 2.07 s and the behavior of the MPC controller was evaluated in three different scenarios:

$$Q\_{\min} \le Q(k+j) \le Q\_{\max} \qquad \forall j = 1, \dots, N \tag{35}$$

• The linear viability condition in Equation (8) is applied, that is, Equation (35), and the following constraint are added to Equation (34):

$$-\frac{\dot{Q}\_{i,\text{min}}}{\dot{Q}\_{i,\text{max}}}Q\_{i}(k+j) + \dot{Q}\_{i}(k+j) + \frac{\dot{Q}\_{i,\text{min}}}{\dot{Q}\_{i,\text{max}}}Q\_{i,\text{max}} \le 0 \qquad \forall j = 1,\dots,N, \quad \forall i = 1,\dots,6; \tag{36}$$

where the subscript *i* denotes the *i*th element of the vector.

The viability constraints are imposed throughout the whole control horizon for the sake of clarity and because the small increment in the computational time (due to the larger number of linear constraints) does not represent a significant issue in the presented experimental tests. However, in the case of linear MPC, imposing such constraints only on the last time instant of the horizon would be enough to ensure feasibility of the control problem and would not give significant differences in terms of control performance. The imposition of the constraint throughout the whole horizon permits to obtain a more reliable prediction of the future trajectory. This can be helpful especially in the case the prediction is utilized in linearized methods such as [15], as a non-reliable prediction might worsen the performance of the method.

Notice that the trajectory is devised in such a way that the desired motions of the first and the sixth joints exceed the maximal viable set, as shown in Figure 7a,b. The figures also show the phase plot of the position/velocity variables computed by the MPC in the three above-mentioned cases. In the case no position bounds are implemented (dashed green line), the joint position obviously does not satisfy the position limit. However, when box position constraints are implemented in the MPC (dash-dotted purple line), infeasibility of the control problem occurs and, since that time, the joint is forced to decelerate with the maximum admissible deceleration, but the position limit cannot be satisfied anyway. This highlights the importance of the viability property. In fact, the implementation of the linear constraints in Equation (8) does ensure the feasibility of the control problem and permits to satisfy the position bound (with a deformation of the original trajectory) (dashed red line). These behaviors are clarified also in Figure 8a,b, where the ideal and the measured joint positions for Joint 1 and Joint 6 are shown.

A second experiment is performed by choosing a total trajectory time *t*end = 2.65 s. In this case, the trajectory drives the first and the third joints in the region between Equations (7) and (8). The behavior of the MPC controller is evaluated in cases where three different viability constraints are implemented:

$$\left(Q\_i(k+j), Q\_i(k+j)\right) \in \mathcal{P}\_{h=3}^i \qquad \forall j = 1, \dots, N, \quad \forall i = 1, \dots, 6; \tag{37}$$

where the superscript *i* denotes the viable polyhedral of the *i*th joint.

• The viable polyhedron obtained by solving Equation (23) with *h* = 4, that is, the following constraints are added to Equation (34):

$$\left(Q\_i(k+j), \dot{Q}\_i(k+j)\right) \in \mathcal{P}\_{h=4}^j \qquad \forall j = 1, \dots, N, \quad \forall i = 1, \dots, 6. \tag{38}$$

Figure 9a,b shows the state trajectory for the first joint in the three different cases. Notice that the position limit is respected and the control problem remains feasible for all cases. However, the reduction of the admissible region determined by Equation (8) gives rise to a significant modification of the original trajectory, although the desired states are always potentially realizable by the robot (as the trajectory does not exceed Equation (7)). An improvement is obtained when the polyhedron with *h* = 3 is used.

Notice that the small velocity bumps visible in the figures are due to predictive nature of the controller. In fact, as the control scheme is based on the tracking of a position reference along the predictive horizon, the controller slightly increases the velocity when it realizes the viability constraint will be activated and a deformation of the task will arise (as the state will be forced to follow the constraint). By giving a small acceleration before the activation of the constraint, the controller minimizes the future deviation with respect to the given position reference.

Finally, the use of *h* = 4 permits to obtain a viable polyhedron that is large enough to enclose the whole desired trajectory. These behaviors are clarified also in Figure 10a,b, where the ideal and the measured joint positions for Joint 1 and Joint 6 are shown.

**Figure 7.** State trajectory in the (*q*, *q*˙)-plane for (**a**) Joint 1 and (**b**) Joint 3 for *Ttot* = 2.07 s: Reference trajectory (gray solid line); MPC without position bounds (green dashed line); MPC with box position bounds (purple dash-dot line); MPC with linear viability inequality (Equation (8)) (red dashed line); maximal viable set (black dotted line); and linear viability inequality ((Equation (8)) (blue dotted line).

**Figure 8.** Position of (**a**) Joint 1; and (**b**) Joint 3 for *T*tot = 2.07 s: Reference trajectory (gray solid line). Output of MPC given as reference signal to the low-level controller: MPC without position bounds (green dotted line); MPC with box position bounds (purple dotted line); and MPC with linear viability inequality (Equation (8)) (red dotted line). Measured position: MPC without position bounds (green dashed line); MPC with box position bounds (purple dash-dot line); and MPC with linear viability inequality (Equation (8)) (red dashed line).

**Figure 9.** State trajectory in the (*q*, *q*˙)-plane for (**a**) Joint 1 and (**b**) Joint 3 for *Ttot* = 2.65 s: Reference trajectory (gray solid line); MPC with linear viability inequality (Equation (8)) (red dashed line); MPC implementing the viable polyhedron with *h* = 3 (purple dash-dot line); and MPC implementing the viable polyhedron with *h* = 4 (dashed green). Viable sets: Maximal (black dotted line); Linear viability inequality (Equation (8)) (blue dotted line); Polyhedron with *h* = 3 (light blue dotted line); and Polyhedron with *h* = 4 (gray dotted line).

**Figure 10.** Position of (**a**) Joint 1 and (**b**) Joint 3 for *T*tot = 2.65 s. Reference trajectory (gray solid line). Output of MPC given as reference signal to the low-level controller: MPC with linear viability inequality (Equation (8)) (red dotted line); MPC implementing the viable polyhedron with *h* = 3 (purple dotted line); and MPC implementing the viable polyhedron with *h* = 4 (green dotted line). Measured position: MPC with linear viability inequality (Equation (8)) (red dashed line); MPC implementing the viable polyhedron with *h* = 3 (purple dash-dot line); and MPC implementing the viable polyhedron with *h* = 4 (green dashed line).

#### **6. Conclusions**

In this paper, we have proposed a method to compute the maximal viable polyhedron for a robot manipulator with bounded joint positions, velocities, and accelerations. In the proposed approach, each joint is considered separately and two optimization problems are devised: one for the upper bounds and one for the lower ones. Given its number of sides, the resulting polyhedron maximizes the area of the admissible position/velocity region. Moreover, the set is proven to be viable and convex. In this way, it can be easily implemented as linear constraints in optimization-based control methods, such as linear MPC. Numerical results demonstrate that the percentage of area covered by the polyhedrons increases as the number of sides grows. This allows the controller to exploit a larger admissible state region and, thus, gives it a broader margin of maneuver in the case the desired motion drives the robot in proximity of the viability boundary. In this case, a better task following can be achieved, as demonstrated by means of experimental results on a six-degrees-of-freedom manipulator.

**Author Contributions:** M.F. devised the method and wrote the paper; M.B. and N.P. conceived and performed the experiments; and A.V. supervised the whole research work.

**Funding:** The research leading to these results was partially funded by the European Union H2020-CS2-738039-Eureca: Enhanced Human Robot cooperation in Cabin Assembly tasks.

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

#### **References**

© 2018 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Motion Investigation of a Snake Robot with Different Scale Geometry and Coefficient of Friction**

#### **Naim Md Lutful Huq, Md Raisuddin Khan \*, Amir Akramin Shafie, Md Masum Billah and Syed Masrur Ahmmad**

Department of Mechatronics Engineering, International Islamic University Malaysia, Selangor 53100, Malaysia<EMAIL_ADDRESS>(N.M.L.H<EMAIL_ADDRESS>(A.A.S<EMAIL_ADDRESS>(M.M.B<EMAIL_ADDRESS>(S.M.A.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +6-03-6196-4575

Received: 1 March 2018; Accepted: 19 April 2018; Published: 28 April 2018

**Abstract:** Most snakes in nature have scales at their ventral sides. The anisotropic frictional coefficient of the ventral side of the snakes, as well as snake robots, is considered to be responsible for their serpentine kind of locomotion. However, little work has been done on snake scales so far to make any guidelines for designing snake robots. This paper presents an experimental investigation on the effects of artificial scale geometry on the motion of snake robots that move in a serpentine manner. The motion of a snake robot equipped with artificial scales with different geometries was recorded using a Kinect camera under different speeds of the actuating motors attached to the links of the robot. The results of the investigation showed that the portion of the scales along the central line of the robot did not contributed to the locomotion of the robot, rather, it is the parts of the scales along the lateral edges of the robot that contributed to the motion. It was also found that the lower frictional ratio at low slithering speeds made the snake robot motion unpredictable. The scales with ridges along the direction of the snake body gave better and more stable motion. However, to get the peg effect, the scales needed to have a very high lateral to forward friction ratio, otherwise, significant side slipping occurred, resulting in unpredictable motion.

**Keywords:** snake robot; snake scale; scale geometry; friction ratio; serpentine motion

#### **1. Introduction**

Snakes possess a unique feature of locomotion that no other creature has. Although legs and wheels give very effective and efficient motion, this creature has the unique ability to move through terrains that are almost impossible to travel through by limbs or wheels. Thus, the ability to move through almost all kinds of surfaces puts snakes in a superior position in terms of locomotion. Snakes have four basic types of locomotion: concertina, serpentine, sidewinding, and rectilinear. Among these, the serpentine locomotion is a faster and a very efficient one for normal conditions. Thus, this motion can be proven exceptionally beneficial in snake robot applications such as search and rescue [1,2], firefighting [3,4], assisting in surgery [5,6], surveillance [7,8] exploring unknown territory, and so forth.

The mechanics of serpentine locomotion were first discovered and explained by J. Gray [9] in 1946. According to Gray, the serpentine locomotion of snakes is due to the lateral friction of the scales of a snake that work as supporting pegs. After the ground breaking discovery of J. Gray, Hu et al. [10] carried out several investigations on real snakes and established the theory experimentally. Some more works on biological snakes that were not focused on the analysis or investigation of the forces applied by the snakes took place. Jayne [11,12] was more focused on the muscular activities during the different types of locomotion with extensive work on serpentine and sidewinding locomotion. Jayne provided extraordinary information regarding the optimization of the actuation method and the energy for the future developments of robots. Meanwhile Miller [13] simulated the motion dynamics

of snake locomotion. Later on, Jayne and Davis [14] developed a relationship between the tunnel width and the speed, or, in other words, the movable space available and the speed of the snake using the kinematics of the concertina locomotion of snake. Marvi and his group studied the friction augmentation by the snakes in concertina locomotion [15].

The archive of information on snake robots dates back to the year 1993 when S. Hirose developed the idea of the locomotion of a snake robot [16], encouraging researchers in working consistently to improve and optimize their artificial snakes. For instance, Wang, Osborne, and Alben [17] and Jing and Alben [18] optimized locomotion on an inclined surface. Instead of real snakes, Kyriakopoulos, Migadis, and Sarrigeorgidis [19] worked on the kinematics of snake robots and provided a design and motion planning in 1999. Prautsch and Mita developed a theory for the dynamic position control [20] for snake robots. Khan et al. [21] first developed a scale based snake trying to mimic the serpentine motion. Later on, they put forward their previous work by investigating the snake robot moving with the serpentine method of locomotion [22]. Apart from that, Marvi et al. [23] presented a snake inspired robot segment based on concertina locomotion in 2011. It was capable of changing the angle of attack of its scales with the change in the slope angle while climbing uphill. They also worked on sidewinding with minimal slips for both the snake and the snake robot on sandy slopes [24]. Varesis, Diamantopoulos, and Tzes [25] conducted an experiment on robots with serpentine locomotion moving through an inclined surface [23]. However, none of the works on artificial snake scales considered the geometry of the scales to investigate the motion of the snake robots. As a continuation of the above efforts, this paper presents the effects of the scale geometry and the friction factor on the locomotion of a snake robot.

This paper is organized as follows. The experimental setup and data acquisition system are described in Section 2, while the effect of the surface properties on motion and the effect of the surface on the directional stability are presented in Sections 3 and 4, respectively. Section 5 draws the conclusions of the paper.

#### **2. The Experimental SETUP and Data Acquisition**

A snake robot made of nine links, eight servo-motors, and a microcontroller, as shown in Figure 1, was used in this research. The microcontroller was programmed for controlling the motion of the servo-motors and propelling the snake in a serpentine motion. The snake's scales were attached at the ventral side of the snake robot, while markers were attached on the dorsal side to keep track of the motion of the robot.

**Figure 1.** The snake robot with markers on the dorsal side.

Four different types of scales were used to investigate the effect of scale geometry and the coefficient of friction on the motion of the snake robot. Photographs and the geometry of the scales are shown in Figure 2a–d. On an average, the span of the scales was 20 mm, except for the scale in Figure 2d. The scales shown in Figure 2a were fabricated from Acrylonitrile Butadiene Styrene (ABS) using a 3D printer while the scales shown in Figure 2b,c were made of ordinary plastic shells and externally threaded cylinder (Cylinder\_1), splitting the cylinders into two halves. In Figure 2d, a scale

made from an externally threaded cylinder with a lower friction factor than Cylinder\_1 (we call it Cylinder\_2) extending between the left and right edge of a link is shown.

**Figure 2.** The photographs of the snake scales attached to the links of the snake robot and orthographic views of the scales (**a**) Designed scale; (**b**) Unthreaded scale (half-shell); (**c**) Threaded scale (Cylinder\_1); (**d**) Threaded scale (Cylinder\_2).

As the goal of this research is to investigate the effect of scale geometry on the motion of snake robots, the robot was equipped with one type of scale at a time. One pair of scales were attached at the lateral edges of each link of the robot. The snake robot was run on three different floor surfaces, namely, Surface 1 (Cloth), Surface 2 (artificial lather 1), and Surface 3 (artificial leather 2). The static friction coefficients between the snake scales and the surfaces are presented in Table 1. One of the main features of the snake scales was the anisotropic friction coefficients along the tangential and normal directions. Frictional anisotropy were achieved through the design of the snake scales of the snake robot.

**Table 1.** The static friction coefficients between the scales and different surfaces.

Figure 3 illustrates the robot orientation and the positions of the scales L1 through R4. Where 'L' denotes the left side and 'R' denotes the right side looking from the top of the robot. A Microsoft Kinect XBOX was attached on a frame above the plane of the motion to capture the motion and the shape of the robot during its motion. A photograph of the experimental setup is shown in Figure 4.

**Figure 3.** The snake's shape and scale placement (L and R stands for left and Right respectively).

**Figure 4.** The experimental setup for the snake's motion capture.

#### **3. The Effect of the Surface Properties on Motion**

The otion data of the snake robot are presented in Table 2.

where, Cycle/s = the fraction of a complete cycle completed by the robot; Dis/Cy = the distance traveled by the robot in one cycle; *V* = the average velocity of the robot;between the initial and final orientations of the robot; High, medium, slow = the cycle speed of the robot.

*Robotics* **2018**, *7*, 18

Table 2 compares four motion elements like the number of cycles per second (seven cycles were taken for each observation), the distance traveled by the robot per cycle, the average velocity of the robot, and the deviation from the initial axis or orientation of the robot after each cycle corresponding to three different motor speeds, namely, high speed, medium speed, and low speed. The three types of motions (high, medium, and low) of the robot were the cycle speed of the sine wave those were created by controlling the speed of the 8 actuating motors. The symbol 'X' in the first row of this table means no data, that is, the robot was unable to move without damaging itself on surface 1 with the designed scale as the scale got stuck in the surface. The velocity of the robot should depend on the frequency of the sine wave of the snake robot. This is evident when the velocities of the robot under three types of motion are compared under the three types of motion. The highest velocity (2.918 cm/s) was achieved in the case of the cylinder\_1 on surface 1 corresponding to the high speed of the motors. This surface combination also produces the highest velocity corresponding to the medium and low speed of the motors. The angular deviation of the robot *θ* from its initial orientation, in the case of the high and medium speeds, are 5.638◦ and 8.607◦, respectively, which are quite high compared to the 0.533◦ that happened in the case of low speed. The friction ratio for this surface combination is lower than that of the cylinder\_1 on surface 2, however, the coefficient of friction along the normal direction for cylinder\_1 on surface 1 is the highest (*μ* = 0.399). Thus, it may be concluded that the higher value of the coefficient of friction along the normal direction plays a significant role in increasing the speed of the robot.

Among the four types of scales and the three types of surfaces, the friction coefficient of surface 2 gives the maximum possible friction ratios between the tangential and normal friction coefficient (Table 1), while the threaded cylinder\_1 on surface 2 gives the highest friction ratio among them all. The comparison of the robot velocities of the designed scale on surface 2 and that of cylinder\_1 on surface 2 showed very similar results. The robot velocities achieved in these cases are almost half that of cylinder\_1 on surface 1 for all types of motor speed. However, the angular deviations are significantly low except for the low speed of the motors. The half shell on surface 3 has the highest movement with better stability (lowest angular deviation) than others, while the smooth half shell gives the minimum distance per cycle. It is noticed that the half shell had the highest friction ratio on surface 3. The last piece of data is a very important one, this data is for the threaded cylinder\_2 on surface 2 (the one with the highest friction ratio) where the scales are placed at the mid surface (centerline) of the links instead of at the edges of the links. It shows that it does not move even a centimeter, yet, the same scale on the same surface gave the best output among all other combinations. Thus, it can be concluded that snake scales need to be placed at the edges of the links for the better motion of the robot, and it is unnecessary to place the snake scales along the centerline of the links.

#### **4. The Effect of the Surface on Directional Stability**

The poses of the robot after each cycle of motion are presented in Figure 5. Each line in these figures is the connecting line between the head and tail of the robot after each cycle of actuation. Thus, these lines are presenting the change in position as well as change in the orientation with respect to the *x*-*y* axis. The initial orientation of the robot was parallel to the *x* axis.

It is observed in Figure 5a–l that the surfaces and the scales have profound effects on the directional stability and motion for the movement of the robot.

Due to practical reason it was not possible to run the robot on surface 1 using the designed scale, thus, the data is absent. However, some important phenomena can be observed from the rest of the figures. All the scales that have threads on it (other than Figure 5c–e,l) are superior in terms of stability in the motion direction, as well as the distance covered. On top of that, when they move on the surface 1 (cloth) (Figure 5f,i) which has an anisotropic friction property within itself, cover higher distance. On the other hand, the half shell (Figure 5c–e) did not have much of a promising outcome in terms of the directional stability and velocity. On surface 2, the motion direction is better maintained than that on surface 3.

**Figure 5.** *Cont*.

**Figure 5.** *Cont*.

**Figure 5.** The poses of the snake robot on the different surfaces at the end of each cycle of actuation.

All the orientations of the snake robot after the last cycle of actuation found in Figure 5a–l are presented in Figure 6. A very important finding that is observed in the last set of data in Figure 5l is shown circled in red in Figure 6. This is the situation when the scales were placed in the middle of the snake body leaving the edges of the link free and floating. There is no significant movement as well as no control of direction in this case. In fact, the robot moved in the backward direction instead of moving forward, however, the same threaded cylinder 2 with the scale at the outer edges of the links gave better robot motion.

**Figure 6.** The comparison of the motion paths of the snake head on different surfaces.

#### **5. Conclusions**

This work makes an effort to discover the relationship between the different variables of the snake scales, like friction ratio, scale geometry, and the span of the scales beneath a snake robot that undergoes serpentine motion. The comparison of the motion shows that scales with a higher friction ratio set at the outer edge of the snake robot helps to generate better and more stable motion. On the other hand, scales only at the center of the belly, irrespective of the friction ratio, does not help the robot to generate any good motion or any forward motion at all. Scales have a ridge like texture on it along the direction of the snake body, which, when placed at the outer edge of the belly, helps generate stable and faster motion. It is predictable from the above experiments that the higher friction ratio will help achieve higher and better motion. Motions with the actuating motors at lower cycle speeds with a scale-surface pair having a lower friction ratio are very unpredictable in nature.

**Author Contributions:** Writing –Original Draft Preparation, Naim Md. Lutful Huq; Supervision, Funding Acquisition & Project administration, Md. Raisuddin Khan and Amir Akramin Shafie; Software, Md. Masum Billah; Methodology, Syed Masrur Ahmmad.

**Funding:** This research was funded by Ministry of Higher Education (MOHE) grant number FRGS13-072-0313.

**Acknowledgments:** The authors profoundly acknowledge the Ministry of Higher Education (MOE), Malaysia for funding this research through the Fundamental Research Grant Scheme (FRGS).

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

#### **References**

© 2018 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Fully Mechatronical Design of an HIL System for Floating Devices**

#### **Hermes Giberti 1,\*, Francesco La Mura 2,\*, Gabriele Resmini 2 and Marco Parmeggiani 2**

Received: 22 May 2018; Accepted: 11 July 2018; Published: 20 July 2018

**Abstract:** Recent simulation developments in Computational Fluid Dynamics (CFD) have widely increased the knowledge of fluid–structure interaction. This has been particularly effective in the research field of floating bodies such as offshore wind turbines and sailboats, where air and sea are involved. Nevertheless, the models used in the CFD analysis require several experimental parameters in order to be completely calibrated and capable of accurately predicting the physical behaviour of the simulated system. To make up for the lack of experimental data, usually wind tunnel and ocean basin tests are carried out. This paper presents a fully mechatronical design of an Hardware In the Loop (HIL) system capable of simulating the effects of the sea on a physical scaled model positioned in a wind tunnel. This system allows one to obtain all the required information to characterize a model subject, and at the same time to assess the effects of the interaction between wind and sea waves. The focus of this work is on a complete overview of the procedural steps to be followed in order to reach a predefined performance.

**Keywords:** parallel kinematic machine; kinematic optimisation; hardware in the loop; mechatronic design

#### **1. Introduction**

Recent simulation developments in Computational Fluid Dynamics (CFD) have widely increased the knowledge of fluid–structure interaction. This has been particularly effective in the research field of floating bodies such as offshore wind turbines and sailboats, where air and sea are involved. Scale model testing of floating structures is a common practice in this research field and help in the development of new concepts and technological solutions, driving the design choices and rapidly answering the scientific questions. Moreover, this study approach allows one to evaluate related problems such as the definition of the control algorithm structures and quantification of costs and benefits.

Nevertheless, when the effects of wind and wave loads become comparable, the validation of the test models is affected by a scaling issue called Froude–Reynolds conflict, related to the ability to reproduce the effects of wind and wave at the same time in a correct manner. For this reason, hybrid tests both in ocean basins [1] and wind tunnels [2] are the most effective approach used to overcome scaling constraints and exploit separately wind and wave generators.

In this paper, the hardware/software setup developed for hybrid tests in a wind tunnel is presented and fully described. The aim of this work is to provide a complete overview of the design procedures and methodologies required to develop a hardware in the loop device capable of simulating the effects of the sea on a physical scaled model positioned in a wind tunnel. In order to achieve the required performance and in compliance with the dimensional constraints, a fully mechatronical design approach has been used. Moreover, this system is particularly challenging in terms of dynamic response. For these reasons, it can be taken as an example of integrated design not only for this specific application, but also for a generic simulator. In fact, the design scheme followed and the numeric tools to develop the system sizing and select the commercially available components are of general interest.

Figure 1 summarises the design procedure of the HIL device presented in this paper and called Hexafloat. The latter consists of a six-degree-of-freedom (DOF) parallel manipulator custom designed to move a scaled wind turbine (e.g., 1/75 of the 10 MW DTU reference wind turbine [3]) in a wind tunnel. The manipulator end effector is equipped with a six-degrees of freedom weight scale used to measure the constraint reactions between the robot and the scaled wind turbine (i.e., aerodynamic load effects). A real-time mathematical model generates the reference commands to be followed by the robot according to the dynamic of the floating substructure simulated by the system and the measures acquired from the weight scale (i.e., real-time hydro/structure computations). The design procedure, as shown in Figure 1, is made up of four main steps that for descriptive simplicity are presented in sequence in this paper but should be considered as a cyclic procedure that ends asymptotically towards the result. The first step relates to the kinematic topological definition and its optimisation as a function of the workspace dimension and the dimensional constraints. The second and the third steps are related to the sizing procedures necessary to ensure the static and dynamic performance of the device taking into account both the actuation system and structural components. The last step deals with the definition of the hardware and software control issues.

**Figure 1.** General scheme of the HexaFloat design process.

#### **2. Geometric and Kinematic Design**

The starting point of every project is the definition of the operating conditions. Figure 2 summarises the dimensional constraints and shows the position of the desired working volume within the wind tunnel. This position has been chosen in order to optimize the layout of the aerodynamic tests using a scaled turbine with respect to the wind tunnel performance. No commercial solutions suit the requirements due to the restricted space available to install the device so therefore a custom solution must be realized out of necessity.

**Figure 2.** Dimensional constraints: (**a**) Workspace constraints; (**b**) Wind tunnel frontal view and dimensions of the space dedicated to the machine.

Among the possible six-DoF manipulator architectures, a parallel kinematic machine (PKM) has been preferred rather than the serial one, due to its high stiffness and high dynamic performance even if this entails an accurate design of the joints in order to have an adequate working space and the absence of backlash. In the first part of the design procedure, the main focus is on the choice of the kinematic topology and the geometric dimensions [2,4]. The most common and widespread six-DOF PKM architectures among which to choose the right topology for this application can be summed up in three families:

The most promising solution within the three families is the six-PUS configuration because the actuation system lays on the floor grounded; moreover, the vertical bulk is limited as opposed to the horizontal one. This topology has been selected in order to guarantee high performance in terms of dynamic response and structural properties with the view to reduce the size vertically. Within this category, the two manipulators shown in Figure 3 have been chosen and evaluated [5]. The first one is called Hexaglide, having parallel linear guides and therefore is characterized by symmetry with respect to the longitudinal median plane. The second is called Hexaslide and is characterized by a radial symmetry with respect to the vertical axis due to the radial distribution of the linear guides in the plane. This characteristic leads to a better isotropy of the workspace.

#### *2.1. Kinetostatic Optimization*

An optimization process [5] is required in order to synthesize the geometric parameters of the two chosen robot architectures and evaluate which is the best solution. To properly set up the process, it is first necessary to identify the geometric parameters characterizing the two architectures and thereafter a function that mathematically describes the goal to be achieved. Finally, a set of physical constraints affecting the system have to be identified and defined in mathematical form in order to create the numeric boundaries for the solution. When choosing the number of parameters to be optimized, a trade-off has to be found: if the number of parameters used to characterize a machine is increased, the whole process of optimization would benefit, thus allowing a higher freedom. However, the complexity of the machine increases dramatically. A trade-off needs to be found between global performance and structural modularity, which is critical in order to simplify design and optimization steps.

**Figure 3.** Hexaglide robot (**a**) and Hexaslide robot (**b**) simplified geometry scheme.

The Hexaslide architecture can be analysed taking into consideration the following parameters (Figure 4) for the optimization process:

**Figure 4.** Hexaslide geometric parameters, the Cartesian space is represented in meters.

For this architecture, the same kinematic chain is replicated identically for six times.

The Hexaglide architecture, shown in Figure 5, can be analysed considering a higher number of parameters because three different kinematic chains can be identified in the robot configuration. The parameters to be optimized are:

**Figure 5.** Hexaglide geometric parameters, the Cartesian space is represented in meters.

For this architecture, each kinematic chain is replicated twice, one on both sides of the longitudinal median plain of symmetry.

The goal is to optimize robot kinematic performance while ensuring the desired workspace boundaries. In order to define a proper objective function, the desired workspace is discretized in elementary volumes: for each elementary volume, the specific set of parameters under analysis must be allowed to reach its centre point. Moreover, for each centre point analysed, the resulting machine must be able to orientate the TCP with any possible combination of pitch, roll and jaw angles varying inside a predefined range. If that is not the case, the elementary volume is added to the total volume that the current set of parameters is not able to cover. In order to make the computational cost affordable and restrict the checks to a finite amount of angles combinations, the three angular ranges are discretized. The evaluation of the objective-function is performed firstly by fixing the mobile-platform orientation Θ*m* = [*αm*; *βm*; *γm*], and thereafter computing the part of volume that end-effector is unable to reach for these specific orientations, such as:

$$\upsilon\_{nc}(\Theta\_m) = \sum\_{i} \sum\_{j} \sum\_{k} c\_{i,j,k}(\Theta\_m) \Delta v \quad \begin{cases} i = 1, \dots, N\_{\mathbf{x}\_{\mathbf{y}}} \\ j = 1, \dots, N\_{\mathbf{y}\_{\mathbf{y}}} \\ k = 1, \dots, N\_{\mathbf{z}\_{\mathbf{z}}} \end{cases} \tag{1}$$

where the combination of parameters *i*, *j*, *k* unequivocally identify a check-point of the workspace in which *Nx*, *Ny*, *Nz* represent the number of discretization points, respectively, along the *x*, *y*, *z* axes. The variable *ci*,*j*,*k*(Θ*m*) is equal to 1 if the specific parameters prevent the end-effector from reaching the pose defined by *i*, *j*, *k* and Θ*m*, and 0, otherwise. The term Δ*v* is the elementary volume resulting from the discretization of the workspace.

This procedure is repeated for all combinations of pitch, roll and jaw angles and the resulting value for the objective-function *Vnc* is computed as:

$$V\_{nc} = \sum\_{m} \upsilon\_{nc}(\Theta\_{m}).\tag{2}$$

Even if a volume portion is discarded for only one specific orientation, that portion of workspace is excluded for the corresponding geometrical configuration. Using this approach, in the ideal case in which a set of parameters allows the robot to reach any point in the workspace regardless of its orientation, the objective-function would be equal to 0. On the other hand, if the set of parameters does not allow the reaching of a particular point of the workspace, the objective function would assume a value equal to the dimension of the single discrete volume times the number of orientations for which that volume has been discarded.

The kinematic capability of reaching each point with every possible orientation is not in itself enough, but additional constraint definitions are required. The kinematic constraints are defined as follows:

The kinetostatic constraints are defined by the transmission ratio between forces and moments acting on the end-effector and the actuation forces. This transmission ratio for each actuator is computed and the maximum should be lower than a prescribed limit value. As is common for PKMs, this transmission ratio varies considerably within the nominal workspace due to nonlinear kinematics, especially in proximity of the singular configurations. The geometric constraints enforce the respect of the minimum distance between two links and between a link and the mobile platform in order to avoid the problem of self-collision of the component for particular poses.

The solution related to the minimisation of the objective function has been obtained through a single objective genetic algorithm approach [6]. The use of a semi-stochastic search and the evaluation of the performance of different individuals at each iteration makes the process of finding the global minimum of the objective-function to be minimized, easier. The steps characterizing this kind of algorithm are:

• Choice of a sufficiently high number of individuals representing a generation in order to have a significant statistical sample;

The algorithm would keep modifying parameters trying to cover the desired workspace as much as possible and minimizing the defined objective function. In Tables 1 and 2, the lower and upper bounds imposed at the beginning of the optimization and the optimal values obtained are reported.

**Table 1.** Limits and results of the Hexaslide optimization.

**Table 2.** Limits and results of the Hexaglide optimization.

Minimum distance between the links, minimum distance between links and platform, transmission ratio and maximum and minimum strokes of the actuated joint coordinates are mapped throughout the workspace to guide one in the choice of the best architecture. The following conclusions hold respectively for the Hexaslide and Hexaglide.

Hexaslide:

#### Hexaglide

It can thus be concluded that Hexaslide architecture is more compact compared to Hexaglide but reaches higher values of force transmission ratio. The Hexaslide architecture has been chosen due to its lower vertical bulkiness, more compact in plane dimensions, better workspace isotropy and lower position of the workspace centre. All of these features make it possible to install the machine under the wind tunnel floor level, reducing the robot influence on the air flow quality and keeping the turbine farther from the wind tunnel ceiling. Moreover, Hexaslide offers two additional advantages: all the elements that make up the links are the same for all six of the kinematic chains and the radial symmetry simplifies the design process. Hereafter, this machine will be referred to as Hexafloat, while the term Hexaslide will refer to the specific architecture of the robot.

#### *2.2. Kinematics Analysis*

In order to develop the optimisation problem and design the control algorithm of the robot, it is necessary to solve the inverse and forward kinematics. In this section, for the sake of briefness, only the solution of the inverse and forward kinematic problem of the robot architecture chosen is presented.

#### 2.2.1. Inverse Kinematics (IK)

In order to study *Inverse Kinematics*, two different reference frames have been considered, the global one and the local one with origin in the TCP and built into the robot platform. With reference to Figure 6, two different vector closures for each kinematic chain can be set up in order to compute the joint coordinates vector *q*. The first vector closure allows one to determine the absolute position vector **di** of the platform joint *Bi* with respect to a point that is the intersection of the **qi** direction with its orthogonal plane passing through the global origin. Each **si** has a fixed length with different orientations in the *X*0-*Y*0 plane:

$$\mathbf{d}\_{\mathbf{i}} = \mathbf{p} + R\_{TCP} \mathbf{b}\_{\mathbf{i}} - \mathbf{s}\_{\mathbf{i}\prime} \tag{3}$$

where *RTCP* is the rotation matrix defining the platform orientation. This matrix is used to transform the expression of **bi**, constant into a rotating local reference frame, in its equivalent with respect to a zero orientation local frame translating with the platform. Position vector **p** and the offset **si** complete the transformation from the zero orientation local frame to the global one. The second vector closure allows one to determine the position of the i-th slider on the guide:

$$\mathbf{d}\_{\mathbf{i}} = \mathbf{d}\_{\mathbf{i}} - \mathbf{q}\_{\mathbf{i}} \mathbf{\hat{u}}\_{\mathbf{i}}.\tag{4}$$

The magnitude of **li** vector corresponds to the length of the robot leg while **uˆ** represents the guide direction unitary vector. This procedure is identical, except for the orientation of **bi**, **uˆi** and **si**, for all six links of the robot.

**Figure 6.** Vector closures used to solve the inverse kinematics (IK) problem: in blue, the 1st closure while in red the 2nd one.

#### 2.2.2. Forward Kinematics (FK)

The Hexafloat robot will not be equipped with sensors able to directly measure the platform position and orientation, therefore a good quality FK could work as a virtual sensor for the pose of the robot. A high frequency estimation of the actual pose of the TCP would be a valuable feedback tool for the other controller in the HIL setup, but a short calculation time and overall stability, in the case of numerical algorithms, is crucial. Extensive research has been conducted to ascertain analytical methods to solve the FK of PKMs, especially for Gough–Steward configuration, but the pose of the robot has not been expressed in explicit form so far.

Considering numerical methods, the Newton–Raphson (NR) algorithm has its numerical stability highly dependent on the accuracy of the initial approximation of the solution vector so a monotonic descent operator can be added obtaining the so-called Modified-Global-Newton–Raphson (MGNR) algorithm, able to estimate the FK solution of six-DoF parallel robots for any initial approximation in the non-singular workspace without divergence. The algorithm requires the definition of a system of six nonlinear equations and the evaluation of a matrix of partial derivatives:

$$
\underline{f}\_i(\underline{X}) = \underline{0}.\tag{5}
$$

$$P\_{ij} = \frac{\partial f\_i}{\partial x\_j}.\tag{6}$$

The evaluation of partial derivatives matrix can be simplified by implementing Jacobian-Free-Monotonic-Descendent (JFMD) algorithm. A first-order Taylor expansion can be used to approximate the partial derivatives matrix in a numerical way. The JFMD method is implemented via the following steps:

$$\mathbf{x}\_{k+1} = \mathbf{x}\_k - \rho\_k [f(\mathbf{x}\_k)]^{-1} f(\mathbf{x}\_k),\tag{7}$$

where *ρk* (0 ≤ *ρk* ≤ 1) is the monotonic descent factor. It starts from 1 and during each iteration is calculated as *ρk* = 2−*m*, where *m* is the number of rechecking times in the corresponding iteration, necessary to obtain the monotonic trend. The approximated Jacobian matrix is evaluated from the first-order Taylor expansion of its partial derivatives and considering a perturbation parameter *η* of 1 × 10−16. The i-th row and j-th column of the approximated Jacobian matrix is evaluated as:

$$J\_{ij} = \frac{f\_i(\mathbf{x} + \eta\_j \epsilon\_j) - f\_i(\mathbf{x})}{\eta\_j},\tag{8}$$

• convergence criterion is defined imposing the error to satisfy the following inequality:

$$\|\|\mathbf{x}\_k - \rho\_k[f(\mathbf{x}\_k)]^{-1}f(\mathbf{x}\_k)\|\|\_2 \le \|\|f(\mathbf{x}\_k)\|\|\_2. \tag{9}$$

• the algorithm stops if convergence is achieved or the maximum number of iterations is reached:

$$\|\|f(\mathbf{x}\_{k+1})\|\|\_{2} \le \delta = 1 \ast 10^{-10} \qquad \text{or} \qquad k \ge k\_{\max} = 20,\tag{10}$$

where *δ* is the required computation tolerance and *kmax* is the given maximum number of iterations.

The first initial approximation must be taken close to the Home Position (463.56 mm from the ground reference frame), where the robot is supposed to be when switched on. For the following times, the initial approximation is chosen as the estimated pose at the previous cycle. The logic of the JFMD algorithm is clarified in the flow chart in Figure 7. Thanks to C language implementation, this routine is able to converge to a solution in relative little iteration, with a mean calculation time of around 50 μs in the actual hardware chosen to control the robot. This performance allows one to use FK as a virtual real-time sensor with sufficient accuracy.

**Figure 7.** Flow chart of the Jacobian Free Monotonic Descendent (JFMD) algorithm.

#### *2.3. Velocity Analysis and Kinetostatic*

In order to determine the relation between the velocity of the TCP and the ones of the joint coordinates, it is necessary to calculate the Jacobian matrix. If all six links are considered together, a compact matrix form coupling the joint velocity vector *q*˙ and the workspace velocity vector *w* can be defined:

$$\begin{bmatrix} \hat{n}\_1^T \hat{n}\_1 & \dots & 0\\ \vdots & \ddots & \vdots\\ 0 & \dots & \hat{n}\_6^T \hat{n}\_6 \end{bmatrix} \dot{q} - \begin{bmatrix} \hat{n}\_1^T & \dots & (b\_1 \times \hat{n}\_1)^T\\ \vdots & \ddots & \vdots\\ \hat{n}\_6^T & \dots & (b\_6 \times \hat{n}\_6)^T \end{bmatrix} w = 0,$$

$$[J\_q] \dot{q} - [J\_{\mathcal{S}^5}]^{-1} w = 0,$$

and with some algebraic steps:

$$\dot{q} = [J]^{-1} w.\tag{11}$$

The solution of the kinetostatic analysis provides the actuation forces *τa* required to bear the external forces *fec* applied to the TCP. Due to the virtual work principle and knowing that the virtual variation of the workspace coordinates is related to the virtual variation of the joint coordinates through the Jacobian matrix [*J*], the actuation forces can be computed as:

$$
\pi\_{\mathfrak{a}} = -[J]^T f\_{\mathfrak{e}\mathfrak{e}}.\tag{12}
$$

The representation of the way in which forces applied to the robot platform are transmitted to the actuators is obtained considering the unitary hypersphere of forces in the workspace [7]. This unitary hypersphere is transformed into a hyper-ellipsoid in the space of actuation forces. This is a common description used in the robotic field in order to characterise the behaviour of the robot in every point of the working space and it states that, if one of the forces applied to the TCP reaches the maximum value of 1, the other components must be nil. An alternative representation is based on the use of a hyper-cube of unitary semi-side that is transformed into a hyper-polyhedron through the Jacobian matrix, as shown in Figure 8. Moreover, it should be noted that the inverse Jacobian matrix can be split into a translational and a rotational part as presented in [8]. Using this approach, it is possible to note that the rotational components of the inverse Jacobian matrix are dimensional and in fact they correspond to a length. In order to let all the elements of the inverse Jacobian matrix be dimensionless, it is useful to divide the rotational components by a scale factor defined as characteristic length *Lc*. In this way, a normalized inverse Jacobian matrix is obtained and is used to obtain the maximum actuation force.

**Figure 8.** Transformation of the workspace forces unitary hypercube into the actuation forces hyperpolyhedron.

#### **3. Actuation Chain Sizing**

A proper sizing of mechanical components and actuating systems requires the knowledge of the most power demanding motion that the machine is required to execute, considering the influence of robot dynamics and carried loads. For the application taken into account, the movement of the end effector is unknown in advance because it depends on the wind effect over the wind turbine and the mathematical model used to define the behaviour of the floating structure. To obtain a scenario of the possible movement of the end effector, it is necessary to consider an approximation of the wave motion of the sea expressed by a combination of six cosinusoidal functions [9,10]. Analysing any combination of a parameter set, which is used to define the sea behaviour, one can obtain a possible end effector motion. The time history of the sea level *ε* is expressed as:

$$
\varepsilon(t) = \rho(t)\cos(\omega t + \varphi(t)),
\tag{13}
$$

where *ρ* is the amplitude, *ω* the pulsation, and *ϕ* the phase shift. Both amplitude and phase shift change very slowly in time so it is reasonable to assume they are constant for the whole duration of the simulation and the time history of each DoF is described as:

$$j(t) = A\_{0,j} + A\_j \cos(2\pi f\_j t + \varphi\_j) \quad \text{with} \quad j = \ge, y, z, a, \beta, \gamma,\tag{14}$$

where *A*0,*j* is an offset parameter taking into account the initial pose of the robot.

The worst operating conditions may be represented in the case where all frequencies *fj* and amplitudes *Aj* assume their maximum value. However, the intrinsically nonlinear kinematic of the robot may invalidate this assumption. In effect, it is not guaranteed that this is the most demanding case in terms of internal loads, motor torques, velocities and accelerations. Furthermore, a simple co-sinusoidal motion of the end-effector is translated into a periodic motion of the joint coordinates where higher order harmonics with respect to the *fj* appear. It is difficult to predict if the energy associated with the higher order harmonics of a specific *fj* is bigger than the energy of the harmonics corresponding to a lower *fj*.

The four characteristic parameters are thus defined in the following ranges:

In order to test any pose in the working volume, a procedure has been designed [11]. Firstly, the workspace has been divided into a finite number of portions. For each of these, an admissible range of motion parameters for each degree of freedom has been chosen. In particular, these ranges have to be discretized in order to obtain a finite number of combinations. Therefore, a set of time histories for the six-DoF end-effector is generated and thereafter the inverse kinematic and dynamic problem solved in order to calculate the motor torques required for each set.

#### *3.1. Multibody Model*

The procedure described above needs a mathematical model to be implemented. In particular, a multibody model of the robot has been developed using the commercial software MSC ADAMS 2016 in order to compute the dynamic and kinematic quantities that allow one to properly size mechanical components as well as the actuating system [12]. At the preliminary stage of the project, the inertial properties of the components, that make up the robot, are unknown and must be estimated. The results obtained by the use of the multibody model and the developed procedure are used to refine these values and update the model to subsequently enhance the results accuracy and converge towards the final result.

The developed multibody model is based on a parametric approach in order to be easily integrated into the motor sizing procedure. The first step for the model formulation is the definition of a set of reference frames located in positions that allow one to easily define the inertial properties, joints, applied forces and measurements points. The reference frames created are: a TCP reference frame; six reference frames called *Bi* in correspondence to the centre of the end-effector joints; six reference frames called *Ai* in correspondence to the centre of the base universal joints; six reference frames called *Pi* located at the origin of the guide axes, as shown in Figure 9.

**Figure 9.** Reference frames configuration, the Cartesian space is represented in meters.

The three main components for each kinematic chain are: platform, links and sliders. The element called platform takes into account the mobile platform, the six spherical joints attached to it, the six axes load cell used to measure the forces and moments exerted by the scaled model of the turbine on the robot. Due to the symmetry, the centre of mass is expected to be located on the *z*-axis of the TCP reference frame and the principal axes of inertia are expected to be parallel to the axes of the TCP. Links can be represented as cylinders whose centre of mass is located in the middle of the link and whose principal axes of inertia are aligned with those of the *Ai* reference frame. The elements called sliders take into account the universal joints, joint supports and sliders of the transmission unit. Since these bodies will be subjected to a purely translational motion during simulations, it is sufficient to characterise them with their total mass. This schematization with the equivalent mass and inertial properties of the assembly is possible because the model described in this section is supposed to be rigid. The mechanical and geometrical limits of each joint are represented by locking one or more relative degrees of freedom and giving a limited range of displacement/rotation among the possible ones.

#### *3.2. Monte Carlo Method*

Due to the impossibility in the definition of the most demanding task for this application, a statistical approach has been developed, based on the Monte-Carlo Method. This novel approach is implemented as follows:

The results collected are used to size the mechanical components and the actuators. In particular, it is possible to obtain the slider velocity and acceleration, forces exerted by the motors on the sliders, axial forces along the guide and the so-called "Load factor" for the actuator sizing procedure. A summary set of values is required in order to easily compare these quantities with the corresponding limits specified by the manufacturers. These sets are computed through the extraction of the maximum values for each simulation stored in *M* array; among these, the highest value is extracted. Therefore, the interval from 0 to the highest extracted value is divided into subranges, and, for each of these, the number of occurrences of the *M* values is assessed. A set of discrete PDFs is obtained and this represents the probability that the maximum value obtained during a simulation is comprised within a certain interval.

#### *3.3. Actuating System Sizing Procedure*

The selection of an actuating system requires one to choose both the electric motor and the gearbox unit. In scientific literature, several procedures to size the motor reducer group are available. In this work, the approach proposed in [11,13] is applied. Independently from the procedure used, the motor reducer sizing is based on the checking of the following three relationships:

• Limit on maximum torque:

$$\max |\mathsf{C}\_{\mathsf{m}}(t)| \prec \mathsf{C}\_{\mathsf{m}, \mathsf{max}} \quad \text{and} \quad d \max |\mathsf{C}\_{\mathsf{m}}(t)| \prec^{\epsilon} \mathsf{C}\_{\mathsf{f}, \mathsf{max}}.$$

• Limits on nominal torque:

$$\mathbb{C}^\*\_{m\mu rms} < \mathbb{C}\_{m\mu nm\nu}$$

• Limit on maximum speed:

$$\max |\omega\_{\boldsymbol{\theta}}(t)| < \tau \omega\_{\boldsymbol{m}, \boldsymbol{\max}} \quad \text{and} \quad \max |\omega\_{\boldsymbol{\theta}}(t)| < \omega\_{\boldsymbol{t}, \max}.$$

*Cm* is the motor torque, *Cm*,*max* the maximum torque the motor is able to generate, and *Ct*,*max* is the maximum torque that the transmission unit is able to bear. In addition, *ωr* represents the resistant speed computed on the load side, *ωm*,*max* is the maximum speed the motor can reach without damaging its mechanical components, *ωt*,*max* is the maximum speed the slider can achieve without damaging mechanical components. Following the [13] approach, from the power balance and the thermal check inequality, it is possible to define the transmission ratio as follows:

$$
\kappa > \beta + \left[ C\_{r,rms}^\*(\frac{\pi\_{rid}}{\sqrt{f\_m}}) - \dot{\omega}\_{r,rms}(\frac{\sqrt{f\_m}}{\pi\_{rid}}) \right]^2,\tag{15}
$$

where the acceleration factor *α* and the load factor *β* are defined as:

$$\alpha = \frac{\mathbb{C}\_{m,nom}^2}{I\_{m}}; \quad \beta = 2 \left[ \left. \dot{\omega}\_{r,rms} \mathbb{C}\_{r,rms}^\* + (\dot{\omega}\_r \mathbb{C}\_r^\*)\_{mcan} \right| . \right. \tag{16}$$

The RMS values of torque and acceleration are defined as:

$$\mathbb{C}^\*\_{r,rms} = \sqrt{\frac{1}{t\_c} \int\_0^{t\_c} [\mathbb{C}^\*\_r(t)]^2 \, dt} \quad \text{and} \quad \dot{\omega}\_{r,rms} = \sqrt{\frac{1}{t\_c} \int\_0^{t\_c} [\dot{\omega}\_r(t)]^2 \, dt}.$$

The acceleration factor *α* is a motor characteristic that can be evaluated directly from the datasheets, whereas the load factor *β* depends on the task performed by the actuations system, and, for the case under study, can be calculated using the multibody model. Following the sizing procedure, firstly, the optimal transmission ratio value *τopt* is evaluated. Then, this value is compared to the transmission ratio purchasable in the market and the closest is chosen taking into consideration that it is within the available range. In particular, greater attention must be paid to evaluate the maximum velocity required during the task and the relative limit transmission value using the formula:

$$\tau\_{lim} = \frac{\max |\omega\_r(t)|}{\omega\_{m,\max}} \max(\tau\_{lim}, \tau\_{\min}) < \tau < \tau\_{\max}.\tag{17}$$

All the results obtained from the multibody simulations refer to the kinematic and dynamic quantities of the manipulator sliders. In order to transform force into torque and linear velocity or acceleration into angular ones, the following equations have to be used:

$$
\dot{\omega}\_{r,i} = \frac{\vec{q}\_i(t)}{\tau\_{TUI}} \quad \mathbb{C}\_{robot,i}(t) = \tau\_{TUI} F\_{robot, \mathcal{X}\_i}(t). \tag{18}
$$

It is worth noticing that to increase the transmission stiffness and achieve the maximum velocity required, every linear axis of the robot is made up of only a single reduction stage. In particular, the motor is directly connected to the screw by means of a rigid joint and the transmission ratio is the screw lead. Due to the high dynamic performance, stiffness and precision required, recirculating ball screws have been chosen. The double slider configuration for each guide has been selected in order to guarantee a uniform load distribution. The manufacturer's instructions to correctly size these devices have been followed, in particular to the total length of the screw, maximum allowed velocities, forces that the guide can bear and the like. All of this is summarized in Table 3. As concerns the electric motor OMRON (Tokyo, Japan), as well as the heat dissipation check, the maximum torque achievable has been evaluated.

**Table 3.** Resume of checks done on guides.

Figures 10 and 11 show the maximum torque distribution and maximum rotation speed obtained during simulations. Once fixed a 99% coverage threshold, the corresponding values are extracted and they are used as guidelines in the components' choice. The most restricting limit is imposed by transmission units: the Rollon ballscrew drive (Vimercate (MB), Italy) chosen satisfies the maximum torque test in 94.7% of the cases with Rayleigh approach and in 89.4% with the other one. Even if the coverage is not 99%, the choice has been made to find a suitable trade-off between performances and size. As for maximum torque analysis, the most restricting limit for rotation speed is imposed by the linear transmission unit, but, since the selected transmission satisfies more than 99% of the cases, this test is passed.

**Figure 10.** Maximum torque distribution along all the simulations.

**Figure 11.** Maximum angular velocity distribution along all the simulations.

#### **4. Mechanical Design and Component Sizing**

A design process is always iterative: an initial model, which satisfies the preliminary parameters required, is firstly created, thereafter each part of this model is refined and improved in order to enhance its dynamic characteristics and reduce mass and costs of the whole system. In this process, the multibody model plays an important role: it is used many times in an iterative way to analyse the dynamic behaviour of the system after every single modification in order to check the mechanical components and actuation system. For the development of a simulator, it is necessary to realize two kinds of structural analyses: the first one is a static analysis for which the focus is to determine the stress and strain of any mechanical components. The second analysis is a dynamic one performed to obtain the modal behaviour of the system.

#### *4.1. Robot Description and Overview*

In order to clarify the results obtained with static and dynamic analysis, a complete overview of the "as built" project is necessary, Figure 12a. The entire system can be divided into two main parts: the Hexafloat manipulator and the auxiliary systems. The first part is made up of "Grounded elements" (all the fixed supporting structure and the power and actuation units), "Joints", "Links" and "Platform", whereas the auxiliary systems include: "Lifting system" (a tool used to manage the rest task of the system), "Energy chain" (device that houses cables and the like) and "External Sensors" (sensors for monitoring the system during the operating tasks).

Hereafter is a more detailed description of every component and technical solution adopted.

**Figure 12.** Components. (**a**) Hexafloat exploded; (**b**) parts on the ground partially exploded view.

#### 4.1.1. Grounded Components

This is the framework and gives stability to all the moving parts, Figure 12b. The linear actuation systems are fixed onto it. Therefore, this assembly should be stiff enough to bear static and dynamic loads. The *Central plate*, made of aluminium, is used to provide the correct orientation and position of the linear guides through calibrated machining. The linear actuators are grouped into two by two and each couple is oriented with an angle of 120 ◦ with respect to the other. This leads to a radial symmetry of the machine. Each couple of linear actuators is fixed to the *Central plate*. To provide further stiffness to the system, avoid undesired tilting and relative displacement between the two couples of guides, three "K" links constrain the *Central plate* with the "C" shaped supports over which the actuation system is mounted. *Aluminium joint holders* are used to couple the *Joints* with the linear actuators. Through the optimization process and the solving of the kinematic closure equations, when the robot is in the *Home Position*, the direction of each *Link* is defined and the *Joints holder* shape and inclination designed. The whole framework will be placed in the wind tunnel using 22 levelling elements to distribute the weight uniformly and maintain the guides on the same plane, without misalignments.

#### 4.1.2. Joints

The PUS kinematic chain on the basis of which the manipulator is realized is made up of an actuated prismatic joint followed by a double revolute joint (universal joint) and a spherical one mounted on the *Platform*. Placing a universal joint and a revolute one so that the rotation axis of the latter passes through the intersection of the universal joint rotation axes, it is possible to obtain the same numbers of DoFs that the spherical joint has. Therefore, the same universal joint used in the lower part of the kinematic chain is mounted on the *Platform*, and the rotational degree of freedom is directly realized within the *Link*. In order to obtain the stiffness, precision and mobility required, a custom joint solution has been realized and they are shown in Figure 13a,b. The main feature of such

*Joint* is the possibility to have a cone shaped motion range with a semi-angular aperture of 45 ◦ around the normal direction, designed and assembled to have zero backlash.

**Figure 13.** Custom joint. (**a**) joint exploded view; (**b**) joint assembly.

Each joint is made up of a couple of *Half shells* connected to an *Inner block* through two roller bearings and two *Support shafts*. These components are tightened by two screws and aligned by two calibrated dowel pins. In this way, the *Inner block* has a relative movement with respect to the shell, the further DoF is provided by another two roller bearings that support the *Inner shaft* pinned to the link. All these inner components are packed through a *Distance ring* and a *Closing ring*. *Half shells* and *Inner block* are made of aluminium, while other components are made of steel. The *Joint* is fixed to the *Joints holder* through four screws and aligned with two calibrated dowel pins.

#### 4.1.3. Links

The *Links*, Figure 14, connect the *mobile platform* to the fixed one. These bear the pay load made up of the mobile platform, *RUAG (RUAG Aviation - Aerodynamic, Schiltwald, Switzerland) load cell*, scaled wind turbine and sensors. The RUAG six-DOF weight scale, model W192–6I, is dimensioned for the following forces and moments: *Fx* = ±1500 N, *Fy* = ±1000 N, *Fz* = ±5000 N, *Mx* = ±500 Nm, *My* = ±1000 Nm and *Mz* = ±600 Nm. Due to the complexity of this device, it is necessary to use a calibration matrix capable of taking into account not only the weight scale deformation but also the working temperature. This calibration matrix, provided by the manufacturer, follows quality procedures in use at the wind tunnel and it is available at the internet address of the Life50+ project. The *Link* is composed of:

**Figure 14.** Link exploded view.

4.1.4. Platform

The *Platform*, Figure 15a, is realised in order to define the right position of the joints and the six-DOF weight scale. It is made up of the following elements:

**Figure 15.** Platform and energy chain. (**a**) platform exploded view; (**b**) energy chain assembly.

#### 4.1.5. Auxiliary System

Two auxiliary systems have been designed: an *Energy chain* (Figure 15b) and a *Lifting system* (Figure 16b). The Hexafloat machine is designed to stay below the floor of the wind tunnel when not in operation. In this configuration, the machine has to overcome a singularity to reach Home Position and a *Lifting system* helps the robot to pass this critical point in the rise and return phases. This system can be schematically modelled as in Figure 16a as an isosceles three-hinged arc, in which hinge A is placed on ground, B on the *Moving carriage* and C at the central part of the *Three points beam*, thus allowing a vertical movement of the *Three points beam* top end (point D). Furthermore, a four-link mechanism is coupled to this system in order to keep it parallel to the ground along the whole stroke of the *Coupling platform*, mounted on top of the lifting mechanism. The *Moving carriage* is driven by the motor through a trapezoidal screw, which converts motion from angular to linear.

Measuring and actuation devices will be mounted on the *Platform* and on-board the scaled wind turbine, and these instruments need to be fed by electric energy and have to be connected to a controller in order to guarantee a real-time data exchange. Therefore, a cable housing system is needed in order to protect cables, maintain a good flexibility while following the *Platform* movements, without interfering with these. The system is made up of the following elements: four "L" plates supports, made of steel; two lateral connecting plates, made of aluminium; two mounting brackets with strain relief mounted, respectively, on the *Platform* and on *Parts on ground*; an intermediate mounting bracket; a mounting sliding bracket and a flexible chain.

**Figure 16.** Lifting system. (**a**) lifting system scheme; (**b**) lifting system assembly.

#### *4.2. Static Analysis*

Through the results obtained from the multibody analysis, one can determine the most critical load acting on the system. In particular, the interest is on the force that must be borne by the links because, through this, the remaining parts of the system are loaded. For safety reasons, an overestimated value of 2500 N is assumed both in traction and compression for every link. The static analysis has been performed using a commercial FEM solution and several models have been studied in order to take into consideration both single elements and their subassembly. Taking into account the final configuration, one can state:

	- **–** Inner Block: a load of 2500 N is split into two equal loads, each one acting on a set of bearings. Maximum stress registered is 85 MPa, which is far below the admissible stress of 250 MPa, Figure 19a.
	- **–** Support shaft: a load of 1250 N is applied for simulating the presence of two Support blocks per Joint. The test is effected by loading one end of the shaft and maintaining the other one fixed, Figure 19b.
	- **–** Inner shaft: a load of 2500 N is applied in the midspan and both ends are pinned in order to simulate the presence of the two bearings, Figure 20a.

inertial load and aerodynamic thrust due to the pay load. These loads are borne very well by the structure both in terms of stress and displacement, Figure 17b.

• Lifting system: this device has to generate a total lift of 150 mm in a time span of about 15 s with a cycloidal motion curve and must also bear the weight of all the moving parts of the robot. This system has been developed in order to help the robot to overcome the singularity configuration during the start and stop procedures. The analysis reveals maximum values occurrence at the initial phase of the rise and FEM analysis ensures that stress remains below the critical values.

**Figure 17.** Movable and fixed platform analysis. (**a**) parts on ground static analysis; (**b**) platform static analysis.

**Figure 18.** Joint analysis. (**a**) joint static analysis; (**b**) joint static analysis, inclined configuration.

**Figure 19.** Inner block and shaft analysis. (**a**) inner block static analysis; (**b**) support shaft static analysis.

**Figure 20.** Inner shaft and link analysis. (**a**) inner shaft static analysis; (**b**) link static analysis.

**Figure 21.** Lower rod analysis. (**a**) lower rod compression analysis; (**b**) lower rod traction analysis.

#### *4.3. Purchased Components' Sizing*

Final choices regarding motors, linear actuators and the energy chain are listed.

#### 4.3.1. Motors

The selected model is the OMRON R88M-K2K030F-BS2, whose main characteristics are collected in Table 4(a). The *Lifting system* motor, required to bear the torque increment given by the low efficiency of the trapezoidal screw adopted, is the OMRON R88M-K20030T-BS2, the main characteristics of which are reported in Table 4(b). Both motor models have an auxiliary brake for safety reasons and an encoder to allow one to have a control of position and velocity. For the 2 KW motor of the robot, the encoder is a quadrature incremental encoder with a maximum resolution of 4,194,304 cnt/rev.

**Table 4.** Motors' characteristics.

The encoder signal is processed in order to obtain a lower resolution of 131,072 cnt/rev. The resulting linear resolution of 3,276,800 cnt/m is obtained considering a lead of 0.04 m/rev for the ball-screw linear axis. For the lifting system motor, the maximum encoder resolution is 131,072 cnt/rev. The modified used resolution is of 13,107.2 cnt/rev, thus obtaining the same linear resolution of the robots axis due to a smaller lead of 0.004 m/rev.

#### 4.3.2. Linear Actuator

The linear actuator chosen is the model TH145 SP4, produced by ROLLON and whose characteristics are reported in Table 5(a). The following customizations are required: screw lead of 40 mm per revolute, to allow the respecting of required performances; two calibrated centring holes on its lower side to correctly assemble the central plate; two calibrated centring holes on the external carriage to allow the correct positioning of the Joint holder.

**Table 5.** Other components' characteristics.

#### 4.3.3. Other Components

The selected energy chain is Igus (Colonia, Germany) Riflex TRE.40 whose main characteristics are reported in Table 5(b). Tapered roller bearings chosen both for the links and joints are the INA FAG 30202-A. These support both axial and radial loads, thus avoiding the usage of more than one bearing. KTR-TOOLFLEX20M motor coupling is chosen. Due to the reduced space available for the *Lifting system*, a more compact model of bearing INA FAG 3000-B-2RS-TVH has been chosen. It is a double row angular contact ball bearing, which can bear both axial and radial loads. KTR Rotex 19/92Sh-A/2.1-*φ*11/2.0-*φ*8 motor coupling has been chosen for the *Lifting system* and it provides the required performance and encumbrance constraints.

#### *4.4. Modal Analysis*

In this section, the approach followed to assess the correctness of the overall design, in terms of dynamic behaviour, is reported. More specifically, the goal is to verify that the first normal modes of the coupled structure coincide with those of the turbine only, in order to minimize the dynamic effects of the robot. Nevertheless, the modal behaviour of a robot strongly depends on the specific pose of the end effector, since the mass and stiffness distribution change according to the pose.

The normal modes and their associated frequencies are computed in order to verify that the frequencies in the first normal mode were high enough. In particular, they must be higher than the well defined frequency range for all the robot poses in to the working volume. The frequency range depends on the dynamic behaviour of the scaled turbine.

To analyse the dynamic response of the manipulator, the entire workspace was discretized, and the trend of the frequency corresponding to the first normal mode over the entire workspace was mapped on specific planes that intersect the robot workspace. The procedure is listed below:

	- **–** Three roll angles: −5◦, 0◦, +5◦,
	- **–** Five pitch angles: −8◦, −4◦, 0◦, +4◦, +8◦,
	- **–** Three yaw angles: −3◦, 0◦, +3◦.

The final result is a set of maps, one for each intersecting plane, which show the trend of the lowest frequencies regardless of the orientation of the robot, as shown in Figure 22. To perform these sets of simulations, a simplified design of the robot was created. Lower simplified joints are tied to the ground applying spherical pins and the simulation provides the first three eigenmode frequencies. Links can be regarded as the main cause for a possible deterioration of the behaviour of the coupled system. The main deformations are concentrated on links and an advanced multibody model has been developed in the "Adams" environment, Figure 23, whose links are made of two extremity elements and an intermediate one characterised by aluminium properties. The results given by the two numerical environments are reported in Table 6.

**Figure 22.** Modal frequencies maps.

**Figure 23.** Numerical modal analysis. (**a**) Inventor 1st mode of vibration; (**b**) Adams 1st mode of vibration.

**Table 6.** Eigenmodes frequencies in Home Position.

A frequency response analysis has been performed on the wind turbine and on the coupled system made by the wind turbine and the Hexaslide. The aim of this simulation is to check the Hexaslide is not influencing the frequencies of the turbine's mode of vibrating, thus making the presence of Hexaslide negligible at least from a mechanical point of view.

A swept sine load of 100 N amplitude between 0.1 Hz and 250 Hz has been applied to the rotor centre in wind (x) direction and the displacements and rotation have been measured. This Force is able to well excite the first and most relevant eigenmodes of the turbine.

The results are reported in Figure 24. The first two wind turbine eigenfrequencies are preserved even in the coupled system,while two peaks appear between 75 Hz and 100 Hz.

It can be concluded that, in the range of 0–60 Hz, a robot's dynamics does not influence the turbine eigenfrequencies.

**Figure 24.** *Cont.*

**Figure 24.** Frequency Response analysis of Turbine and Hexaslide with turbine. (**a**) displacement along *x*ˆ*O* direction; (**b**) rotation around the *y*ˆ*O* axis.

#### **5. Control Architecture and Electronics**

A simplified scheme of the electric panel configuration is reported in Figure 25. The core of the electrical panel is represented by Power PMAC, the property of Delta Tau Data Systems (Chatsworth, CA, USA). There are essentially two zones: one with AC voltage, represented by solid black and blue lines, powered by a 380 V AC line and brings it to the Power PMAC fed at 230 V and to the power circuits of the servo amplifiers fed at 380 V, and the other one with 24 V DC voltage, downstream of the 24 V power supply, which is represented by the solid and dashed green lines. This circuit provides power to the safety relay, limit switches and proximity sensors and to the Beckhoff modules for Ether-CAT communications. Red lines represent the transmission of data between different components.

The modular rack is the most flexible configuration, since it permits the user to choose which CPU card, digital or analogue I/O card, axis interface cards and the like to use in the system. Power PMAC can handle all of the tasks required for machine control, constantly switching back and forth between the different tasks thousands of times per second. On this powerful controller, the main control software of the robot has been designed and implemented. Figure 26 shows the general overview of the software architecture. In particular, it is structured so that a primary and a secondary states machine manages the principal functionality of the system, among which the management of the logic

state (start, stop, homing, jog, and the like) and the functioning of the machine (managing exceptions, motion programs, debug towards the user). The Position-Based-Admittance-Control (PBAC) control scheme used to close the position loop is highlighted in light blue. Further details in this regard are set out in [14].

**Figure 25.** Simplified scheme of an electrical layout.

The lowest level of the control is constituted by the position and velocity servo loops, giving the analogic torque output reference for the seven servo actuator. Each high level control modality developed has been designed to respect real-time performances desired as well as safety, using advanced tools such as position based admittance control, buffering with time-based control, motion look-ahead for smooth blending, fast C written nonlinear FK and demanding algebra, workspace boundaries check

with controlled dynamics on the limits [15], and acceleration saturation with workspace reference tracking check. For more details about control architecture and algorithms, refer to [14,15].

The Human–Machine Interface (HMI) designed allows full control by the user of every significant functioning parameter, simplifying state change and enhancing safety. The HMI communicates with Power PMAC by means of Telnet communication protocol and is completely written in C# language.

**Figure 26.** Complete control scheme.

#### **6. Conclusions**

The experimental verification of the dynamic response of the robot is currently being finalized. The main control features have been tested on a scaled model of the manipulator with optimum results [14,15]. The full scale is now fully operative and its behaviour is under test. For the sake of completeness in Figure 27, the full scale system equipped with accelerometers during the modal analysis campaign is reported. From the first broad results, it appears confirmed that the robot will not interfere, from a dynamic point of view, with the wind turbine scale model [3].

Further experimental modal analysis is staged for the next months, to fully explore the performances not only of the complete mechanical system, but of the whole machine under control influence. First verifications of control performances have already given satisfying results, but further refinement is done every day to an optimized controlled system behaviour in every required operative condition.

This paper shows the design methodology of the HexaFloat system, a six-DoF robot for wind tunnel hybrid testing of floating offshore wind turbines. This setup consists of a parallel kinematic robot, "HexaFloat", designed and developed to test the dynamics of floating offshore wind turbine concepts, selected within LIFES50+ project, at the Politecnico di Milano wind tunnel, through a hybrid methodology which combines HIL, in real-time measurements (i.e., aerodynamic forces on the wind turbine scale model) and computations (i.e., hydrodynamic forces on platform). This represents the complementary test approach, with respect to the one developed at the Sintef Ocean basin. The final test campaign for LIFES50+ project is staged for July 2018. The complete HIL setup is currently under testing for performance verification, disturbance influence analysis, control of final refinements and tuning and controllers' communication optimization.

**Figure 27.** Hexaslide experimental setup.

**Author Contributions:** H.G. and F.L.M. conceived the devices and the sizing tools used; G.R. and M.P. designed the system; F.L.M. developed the control software of the devices; F.L.M., tested the system. H.G. and F.L.M. wrote wrote the paper.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant agreement No. 640741.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**

© 2018 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 (http://creativecommons.org/licenses/by/4.0/).

## *Technical Note* **Nominal Stiffness of GT-2 Rubber-Fiberglass Timing Belts for Dynamic System Modeling and Design**

#### **Bozun Wang 1, Yefei Si 2, Charul Chadha 3, James T. Allison 1 and Albert E. Patterson 1,***∗*

Received: 20 October 2018; Accepted: 19 November 2018; Published: 21 November 2018

**Abstract:** GT-style rubber-fiberglass (RF) timing belts are designed to effectively transfer rotational motion from pulleys to linear motion in robots, small machines, and other important mechatronic systems. One of the characteristics of belts under this type of loading condition is that the length between load and pulleys changes during operation, thereby changing their effective stiffness. It has been shown that the effective stiffness of such a belt is a function of a "nominal stiffness" and the real-time belt section lengths. However, this nominal stiffness is not necessarily constant; it is common to assume linear proportional stiffness, but this often results in system modeling error. This technical note describes a brief study where the nominal stiffness of two lengths (400 *mm* and 760 *mm*) of GT-2 RF timing belt was tested up to breaking point; regression analysis was performed on the results to best model the observed stiffness. The experiments were performed three times, providing a total of six stiffness curves. It was found that cubic regression mod els (*R*2 > 0.999) were the best fit, but that quadratic and linear models still provided acceptable representations of the whole dataset with *R*2 values above 0.940.

**Keywords:** timing belt; belt stiffness; dynamic system modeling; mechatronic systems; 3D printers; robotics

#### **1. Introduction**

Timing belts are a common means of motion transfer between rotating motors/shafts in a machine or mechatronic system. Many small-to-medium sized mechatronic systems such as 3D printers [1], robots [2,3], desktop computer numerical control (CNC) machines [4], and positioners [5] use such belts, typically in the GT-style [6,7]. GT-style belts are specifically designed to effectively translate rotating motion from pulleys into linear motion with minimal deformation, slippage, and backlash. One of the fundamental characteristics of such a motion transfer system is that the length of the belts changes with time, causing time-variant stiffnesses in the belts which must be considered in dynamic system modeling and design. Note that the "stiffness" in the belt is considered only in the tension direction of the belt for this work, resulting in a stiffness that can be described as a single value or function instead of the full stiffness matrix [8,9].

When analyzing and designing any robotic and other mechatronic systems, it is vital that a good dynamic model of the system be developed and used. Since such systems often use some kind of flexible belts for motion transfer, the belt stiffness is a very important parameter in a system model. In cases where the length of the belt is constant (e.g., running between two fixed pulleys), the stiffness

*k* of the belt can be modeled as a spring where *f*(*x*) = *kx*; therefore, the stiffness of the belt is a function of its deflection *x* under load. In effect, this constant-length stiffness of the belt is its "nominal" design stiffness. However, in cases where the belt changes length during use, the effective length of these belt is a function of time and, therefore, its stiffness is also time-variant; this time-variant stiffness is the "effective" or apparent stiffness of the belt at some time *t*. It has been shown that the effective stiffness of the length-changing belts can be directly calculated as a function of the nominal stiffness value, the belt width, and the real-time length of the belt [10] such that:

$$k\_i(t) = \mathcal{C}\_{sp} \frac{b}{L\_i(t)} \tag{1}$$

where *ki* is the effective stiffness as a function of time, *Csp* is the nominal stiffness, *b* is the belt width, and *Li*(*t*) is the length of the belt section at time *t*. For any case where the length remains constant, the effective and nominal stiffnesses are equal since the value of *Li*(*t*) is a constant. Note that the value of *Csp* may be a constant or function of material properties for different belt materials; it cannot be considered a function of time the way that the length of the belt is. The most commonly-used GT-style belt is the GT-2; Figure 1 shows the fundamental geometry and specifications for this type of belt.

**Figure 1.** GT-2 belt (**a**) specifications and (**b**) basic geometry.

Figure 2 shows a common application, where a GT-style belt is used to transfer motion from a stepper motor to drive a linear positioning system. Also shown is a 2D dynamic model representation of such a system (Figure 2b), where the differences in effective stiffness, based on belt length, in the belt sections are clearly evident. The sections *L*1 and *L*2 change in effective stiffness as a function of time, while section *L*3 stays constant during use [11,12] so the effective and nominal stiffnesses are equal.

**Figure 2.** (**a**) simple positioning system that utilizes a GT-type belt to drive the table and (**b**) its representative dynamic model.

The work described in this note explored the nominal stiffness *Csp* and the best way to model it in dynamic systems where belt length is not constant. Several previous studies have assumed that rubber-based timing belts have a linear nominal stiffness [5,11–18]. However, it is vital for designers and engineers working with dynamic systems which use belts for energy transfer to understand the true effects of the belt stiffness [19,20]. Therefore, experimental data was collected and used to derive conclusions on the true stiffness behavior of the GT-2 belts during use. The collected data was subjected to regression analysis to see which type of model best fit, allowing the comparison of models for the same dataset. The information in this study will prove useful, both in choosing *k* stiffness values for dynamic models and for judging expected model error if linear stiffness assumptions are used.

#### **2. Procedure and Results**

Two lengths of new GT-2 belts, 400 mm and 760 mm, were subjected to a simple tensile test until they ruptured. The test apparatus was a custom-built, screw-driven manual desktop test stand set up for tensile testing with 3000 *N* capability and a travel rate of 1/16 *in* (1.6 mm) per screw revolution. The screw drive was rotated at a constant rate of 0.5 revolutions per second (0.8 mm/s), a reading being taken every revolution of the screw or every 1.6 mm. Since the length measurement was based on a count of the threads during travel, the uncertainty in length was too small to quantify; the digital readout for the unit used a load cell with a given uncertainty of 100 gram-force or 0.89 *N*. It was necessary to use this kind of manual tensile testing machine as none of the available standard machines were sensitive enough to measure the force-deflection behavior of these kinds of belts [21]. In addition, the discrete time measurement ensured a reasonably-sized dataset for curve-fitting. This was replicated twice to obtain a set of six different curves, three from each length. The ruptured belts were observed to fail suddenly and to show tearing of the glass fibers inside, as shown in Figure 3. The GT-2 belts used were a composite of neoprene (synthetic rubber) [22] and glass fibers, where the fibers appeared to drive the failure point of the belts.

**Figure 3.** Belt break interface, showing broken fibers.

The collected data, in terms of force-deflection curves, are shown in Figure 4a, while the equivalent stress–strain curves for the tests are shown in Figure 4b. The length of the belts clearly had an effect on the force-deflection curves, but this largely disappeared when the length was accounted for in the stress–strain curves. Note that most of the curves show hyper-elastic behavior, i.e., there is no region in the curve where the stiffness is constant.

**Figure 4.** Belt stiffness curves: (**a**) force-deflection curve and (**b**) stress–strain curve.

As the nominal compliance of the belts was clearly found to be nonlinear, a regression analysis was performed to model the curves and find the level of unexplained variance in these curves. One of the most common polynomial regression models [23,24] used for hyper-elastic materials is the cubic polynomial. The basic model used for this study began with the following polynomial model:

$$
\sigma\_{belt} = A\varepsilon\_{belt}^3 + B\varepsilon\_{belt}^2 + C\varepsilon\_{belt} + D \tag{2}
$$

where a cubic model includes all of the variables, a quadratic model can be generated by setting *A* = 0, and a linear model can be used with *A* = *B* = 0. These curve fits, completed using Microsoft Excel™ (Microsoft Corp, Redmond, WA, USA) are shown in Figure 5a, and the fits for each of the variable and the resulting *R*2 values are shown in the first six cases in Table 1.

After fitting the cubic models to each of the six sets of experimental data, the cubic model, a quadratic model, and a linear model were then fit to the entire set at once, as shown in Figure 5b. A significant drop in the *R*2 value was noted for all of the models fit to the dataset, but differences between the cubic, quadratic, and linear models were observed to be small, as shown in Table 1.

It was observed that the low-strain region of the dataset (Figure 5b,c) conforms better to a linear model when the entire dataset is used. In actual use, it is most likely that the belts will not reach more than 20–30% of the belt breaking strength during normal use [14,25,26], so this is a valid assumption for many systems; this will, of course, need to be determined by the modeler or designer before using a linear belt model. If the low-strain assumption can be used, then the data fit a linear model with a slightly greater *R*2 value than a quadratic model for the entire dataset and is certainly superior to a linear model for the entire dataset. The linear model for this case is shown in the last row of Table 1.