**1. Introduction**

In industrial robotics, stiffness property has grea<sup>t</sup> importance especially when they are used in manufacturing processes. Production quality during milling and drilling operations depends on the stiffness performance of the robot. Thus, many academic studies are devoted to the stiffness analysis and evaluation of industrial robots [1–7].

One reason for the stiffness performance evaluation of a manipulator is to determine the capabilities of manipulators [4,8–11]. There are several evaluation approaches of stiffness properties. One is commonly used in industrial robots that consider the absolute end-effector deflection. In this approach, a specific position or a trajectory is given to the robot, and its positioning error is observed under external wrenches [8,12,13]. The procedures of this kind of evaluation are given in ANSI and ISO standards [4,14]. These methods are suitable to have a general idea on stiffness performance, but they may not distinguish whether forces or moments cause the elastic translational or rotational displacements. Depending on the performance output, a suitable manipulator for the desired task can be selected.

Another reason is to optimally design the manipulator parameters by considering the stiffness performance [5–7,13,15–17]. The optimal design process of a manipulator requires objective functions that include the stiffness performance indices. These indices are mainly used to achieve high rigidity while preserving low inertia and high kinematic performance.

In parallel robots, stiffness modulation can be achieved by introducing redundancy in actuation [18,19] and/or kinematics [20–22]. Stiffness characteristics of the manipulator are regulated by applying suitable redundancy resolution algorithms in control [19–21,23,24]. A common way of redundancy resolution is via null-space control [25] in which a performance index is designed. This performance index is defined as a stiffness performance index when it is required to modulate the stiffness of the manipulator [18–24]. A key aspect of this performance index is that it should be computationally efficient.

Optimal design and control processes require a mathematical stiffness model (stiffness matrix) of the manipulator to use the stiffness performance indices. These indices adopt common matrix operations such as the determinant, norm, or singular value decomposition (SVD) in defining the

performance indices [3–6,8–11,13,14,16,17,26–35]. Each of them investigates a different property of the stiffness matrix. For instance, the determinant is accepted as a scalar indication of stiffness magnitude [3,11,14,16,26,28]. SVD computes the eigenvectors and eigenvalues that show the most/least stiff directions and their magnitudes [3,4,6,8,9,11,13,14,17,26–28,30]. Once the stiffness model is known these approaches are quite handy and they provide detailed stiffness information of the manipulator. However, their general problems are that they are calculated locally.

Isotropy is having the same stiffness performance distribution in all directions. In the literature, eigenvalues of the stiffness matrix are used for relative stiffness resolution evaluation. The evaluation is carried out by calculating the ratio of the maximum to the minimum eigenvalue. The index is named as stiffness condition number, stiffness dexterity, or stiffness isotropy index [10,13,14]. It obtains a directional stiffness resolution compared to the maximum eigenvalue. However, only focusing on maximum and minimum eigenvalues hides the effect of intermediate eigenvalues on the performance. For instance, there exist four possible cases for the distribution of eigenvalues of a 3 *×* 3 stiffness matrix.The first case is the one that has all different eigenvalues. The second case has two equal minimum eigenvalues. The third case has two equal maximum eigenvalues. The fourth case has all equal eigenvalues. For the first three cases, if the minimum and maximum eigenvalues are the same, then, the stiffness isotropy index computation will output the same value. However, it is clear that the isotropy of first, second, and third cases should be different from each other. This problem can only be resolved by taking into account the intermediate eigenvalues in stiffness isotropy computation.

In this study, special attention is directed towards this intermediate eigenvalue problem. In this regard, a new performance index is proposed. This index composes a volumetric isotropy index by considering intermediate eigenvalues. As a case study, an R-CUBE [36] parallel manipulator is considered. A comparison is conducted between the proposed volumetric isotropy index and existing isotropy index. R-CUBE is selected for this analysis since it has relatively trivial kinematics to be in conveniently used in the formulation of its stiffness model [36].

The remaining sections of this paper are organized as follows; in Section 2 stiffness model of the R-CUBE is obtained, and notations are given. Then, in Section 3, currently used performance indices are discussed, and a new performance index is introduced. In Section 4, results of conventional isotropy index and the new index computations are presented. Finally, conclusions are stated in Section 5.

#### **2. Stiffness Model of R-CUBE**

The R-CUBE manipulator is introduced by [36]. The manipulator composed of revolute joints, only. In total, 3 serial chains exist, and each one of them controls their respective translational DoF in Cartesian space. In Figure 1 the kinematic model is illustrated.

In Figure 1a, first refecence frames of the each serial chain are located on *u*(0) *k* orthogonal axes along *k*th axis for *k* = 1, 2, 3. *ij* stands for the *j*th frame in *i*th serial chain as shown in Figure 1b. *u*(0) *k* - *u*(*p*) *k* due the kinematic constraints. *p* is mobile platform frame. Besides, *u*(15) 3 ,*u*(25) 3 ,*u*(35) 3 are always aligned with *u*(35) 1 ,*u*(15) 1 ,*u*(25) 1 vectors, respectively. The forward kinematics is given as:

$$r\_i = S + l\_1 \sin q\_{i1} \text{ for } i = 1, 2, 3 \text{ and } \bar{r} = \begin{bmatrix} r\_1 & r\_2 & r\_3 \end{bmatrix}^T \tag{1}$$

where *S* is a constant distance from 0th frame to *u*(*<sup>i</sup>*0) 3 frame. *r*¯ denotes position vector with respect to the origin.

Commonly used stiffness (elasto-static) modeling methods are classified as Finite Element Method (FEM) [37,38], Matrix Structural Method (MSM) [39,40], and Virtual Joint Method (VJM) [28,41–45]. FEM exhibits the highest accuracy in exchange of computation cost due to its numerical approach. Besides, a meshing operation must be conducted for each pose of the manipulator. VJM and MSM are faster in this evaluation since they construct semi-analytical or analytical models. Thus, stiffness

performance metrics are used via VJM and MSM for evaluation of a manipulator in any of its pose. In this study, VJM approach is adopted.

**Figure 1.** Kinematic sketches of the R-CUBE mechanism: (**a**) The manipulator, (**b**) variables of *i*th serial chain where *i* is one of the serial chains.

Stiffness model of the R-CUBE manipulator is computed by combining the stiffness models of each serial chain. Compliance model of *i*th serial chain and its connection to base and mobile platforms are illustrated in Figure 2. Passive and active joints have 1 degree-of-freedom (DoF) while virtual joints have 6 DoF (3 translations+3rotations). A virtual joint is defined as:

$$H\_{\mathbb{P}}(\bar{\theta}\_{\text{ij}}) = T\_1(\theta\_{\text{ij}}^1)T\_2(\theta\_{\text{ij}}^2)T\_3(\theta\_{\text{ij}}^3)\mathbb{R}\_1(\theta\_{\text{ij}}^4)\mathbb{R}\_2(\theta\_{\text{ij}}^5)\mathbb{R}\_3(\theta\_{\text{ij}}^6) \tag{2}$$

where *Hv* denotes the homogeneous transformation matrix (HTM). *Tk* and *Rk* are pure translational HTM along and pure rotational HTM about *uk*th axis for *k* = 1, 2, 3. ¯ *θij* is a vector that contains virtual joint variables. Superscripts of *θij* represents the element number.

**Figure 2.** Compliant kinematics of the manipulator where AJ is active joint, PJ is passive joint, VJ is virtual joint, MP is mobile Platform, and B is base.

The compliant kinematic model is computed as follows:

$$\begin{aligned} H^{(0,K\_{i1})} &= \mathcal{R}\_3(\boldsymbol{\varrho}\_{i1}) \mathbf{T}\_{1u}(l\_1) \mathbf{H}\_v(\bar{\theta}\_{i1}) \\ H^{(K\_{i1},K\_{i2})} &= \mathcal{R}\_3(\boldsymbol{\varrho}\_{i2}) \mathbf{R}\_1(-\pi/2) \mathcal{R}\_3(\boldsymbol{\varrho}\_{i3}) \mathbf{T}\_1(l\_2) \mathbf{H}\_v(\bar{\theta}\_{i2}) \\ H^{(K\_{i2},K\_{i3})} &= \mathcal{R}\_3(\boldsymbol{\varrho}\_{i4}) \mathbf{T}\_1(l\_3) \mathbf{H}\_v(\bar{\theta}\_{i3}) \\ H^{(K\_{i3},\tilde{s}5)} &= \mathcal{R}\_3(\boldsymbol{\varrho}\_{i5}) \\ H\_{Ki} &= H^{(0,i0)} \boldsymbol{H}^{(0,K\_{i1})} \boldsymbol{H}^{(K\_{i1},K\_{i2})} \boldsymbol{H}^{(K\_{i2},K\_{i3})} \boldsymbol{H}^{(K\_{i3},\tilde{s}5)} \boldsymbol{H}^{(K\_{i3},\tilde{t}5)} \\ H\_{Ki} &= \begin{bmatrix} \mathcal{R}\_{Ki} & \bar{r}\_{Ki} \\ \mathbb{0}^T & 1 \end{bmatrix} \end{aligned} \tag{3}$$

where *HKi* is the compliant transformation matrix of the *i*th serial kinematic chain of R-CUBE, and *HK*1 = *HK*2 = *HK*3 by the assumption of a rigid mobile platform.

*Machines* **2019**, *7*, 44

Active joints are assumed to be rigidly locked to exclude actuation stiffness in the structural stiffness calculations. Then, virtual and passive joint variables are arranged in a matrix form as follows:

$$\boldsymbol{\ddot{q}}\_{pi} = \begin{bmatrix} \boldsymbol{\varrho}\_{i2} & \boldsymbol{\varrho}\_{i3} & \boldsymbol{\varrho}\_{i4} & \boldsymbol{\varrho}\_{i5} \end{bmatrix}\_{4 \times 1}^{T}, \quad \boldsymbol{\ddot{\theta}}\_{i} = \begin{bmatrix} \boldsymbol{\ddot{\theta}}\_{i1}^{T} & \boldsymbol{\ddot{\theta}}\_{i2}^{T} & \boldsymbol{\ddot{\theta}}\_{i3}^{T} \end{bmatrix}\_{18 \times 1}^{T}, \quad \boldsymbol{\ddot{\mathcal{Q}}}\_{i} = \begin{bmatrix} \boldsymbol{\theta}\_{i}^{T} & \boldsymbol{\ddot{\theta}}\_{pi}^{T} \end{bmatrix}\_{22 \times 1}^{T} \tag{4}$$

where *q*¯*p<sup>i</sup>* denotes passive joint variables, *Q* ¯ *i* is generalized coordinates of *i*th serial chain. *Q* ¯ *i* is used to obtain Jacobian matrices as follows:

$$\frac{\partial H\_{Ki}}{\partial Q\_{ik}} = \begin{bmatrix} \frac{\partial R\_{Ki}}{\partial Q\_{ik}} & \frac{\partial r\_{Ki}}{\partial Q\_{ik}} \\ \vec{0}^T & 1 \end{bmatrix} \text{ for } k = 1, 2, \dots, 22 \tag{5}$$

where subscript *k* denotes the *k*th variable of *Q* ¯ *i*.

The Jacobian matrix *JKi* for passive and virtual joints is given as follows:

$$J\_{Ki} = \begin{bmatrix} J\_{Ki1} & J\_{Ki2} & \dots & J\_{Ki22} \end{bmatrix}\_{6 \times 22} \tag{6}$$

where ¯ *JKik* for *k* = 1, 2, ..., 22 denotes the Jacobian column matrices that are related with *k* th variable of *Qi*.

*JKi* can be divided into sub-matrices as *Jθi* and *J pi*. They denote Jacobian matrices of virtual and passive joints. Sub-matrices are presented as follows:

$$\begin{aligned} \mathbf{J}\_{\theta i} &= \begin{bmatrix} \mathbf{\bar{J}}\_{\text{Ki}1} & \mathbf{\bar{J}}\_{\text{Ki}2} & \dots & \mathbf{\bar{J}}\_{\text{Ki}18} \end{bmatrix}\_{6 \times 18} \\ \mathbf{J}\_{pi} &= \begin{bmatrix} \mathbf{J}\_{\text{Ki}19} & \mathbf{J}\_{\text{Ki}20} & \mathbf{J}\_{\text{Ki}21} & \mathbf{J}\_{\text{Ki}22} \end{bmatrix}\_{6 \times 4} \\ \mathbf{J}\_{\text{Ki}} &= \begin{bmatrix} \mathbf{J}\_{\theta i} & \mathbf{J}\_{\text{Fi}} \end{bmatrix}\_{6 \times 22} \end{aligned} \tag{7}$$

Note that, except in a kinematically singular configuration, *Jθi* is always full-rank since the virtual joints have decoupled DoF. Similarly, *J pi* gets rank deficient in kinematic singularities.

Jacobian matrices relate small deflections of joints to task space variables:

$$\begin{aligned} \mathcal{X}\_i &= \bar{f}(\bar{Q}\_i)\_{\prime} \implies \Delta \mathcal{X}\_i = J\_{Ki} \Delta \bar{Q}\_i \\ \Delta \bar{X}\_i &= J\_{\theta i} \Delta \bar{\theta}\_i + J\_{pi} \Delta \bar{q}\_{pi} \end{aligned} \tag{8}$$

where Δ*X* ¯ *i* is 6 *×* 1 column matrix of translational and rotational compliant deflections in task space that are calculated for the *i*th serial chain. Δ operator denotes the change between initial and final states. In the initial state, there is no applied wrench. In the final state, an external wrench is applied to the manipulator.

The relation between the external force/torque applied at the end-effector and joint space force/torque is provided via a property of Jacobian matrix that is given as follows:

$$\mathbf{F}\_{Ki} = \mathbf{J}\_{Ki}^T \mathbf{F}\_{ext} \tag{9}$$

where [*F*¯*Ki*]<sup>22</sup>*×*<sup>1</sup> is the joint space force/torque vector. [*F*¯*ext*]<sup>6</sup>*×*<sup>1</sup> is the external wrench. *<sup>F</sup>*¯*Ki* is divided into sub-components and the force/torque vector on each joint are found as follows:

$$\begin{aligned} \bar{F}\_{Ki} &= \begin{bmatrix} F\_{\theta i}^T & F\_{pi}^T \end{bmatrix}^T\\ \begin{bmatrix} F\_{\theta i}^T & F\_{pi}^T \end{bmatrix}^T &= \begin{bmatrix} J\_{\theta i} & \mathbf{0} \end{bmatrix}^T \bar{F}\_{ext} + \begin{bmatrix} \mathbf{0} & J\_{pi} \end{bmatrix}^T \bar{F}\_{ext} \end{aligned} \tag{10}$$

where *F* ¯ *θi* and *F* ¯ *pi* are force/torque vectors of virtual and passive joints. *Machines* **2019**, *7*, 44

> *<sup>F</sup>*¯*Ki* is a function of the stiffness matrix and deflections that are defined in joints space.

$$F\_{\rm Ki} = \text{diag}(\mathbf{K}\_{\theta i}, \mathbf{K}\_{\rm pi}) \Delta Q\_i \text{ and } \mathbf{K}\_{\theta i} = \text{diag}(\mathbf{K}\_{\theta i1}, \mathbf{K}\_{\theta i2}, \mathbf{K}\_{\theta i3}) \tag{11}$$

where *Kθik* denotes stiffness matrix of *k*th link as expressed in [46] for a beam. Parameters of this stiffness matrix depends on the link geometry and material property. *Kθi* denotes the stiffness matrix of *i*th serial chain. *<sup>K</sup>p<sup>i</sup>* denotes stiffness of passive joints. A relation is expressed in Cartesian space as follows:

$$\begin{aligned} \boldsymbol{J}\_{\bar{K}i}^{\top} \bar{\mathbf{F}}\_{\text{ext}} &= \text{diag}(\mathbf{K}\_{\theta i}, \mathbf{K}\_{pi}) \Delta \bar{\mathbf{Q}}\_{i} \text{ and } \boldsymbol{J}\_{\bar{K}i}^{-1} \Delta \bar{\boldsymbol{X}}\_{i} = \Delta \bar{\mathbf{Q}}\_{i} \\ \Rightarrow \bar{\mathbf{F}}\_{\text{ext}} &= (\boldsymbol{J}\_{\theta i}^{-T} \mathbf{K}\_{\theta i} \boldsymbol{J}\_{\theta i}^{-1} + \boldsymbol{J}\_{pi}^{-T} \mathbf{K}\_{pi} \boldsymbol{J}\_{pi}^{-1}) \Delta \bar{\boldsymbol{X}}\_{i} \end{aligned} \tag{12}$$

Passive joints do not generate reaction torques about their rotation axes. Hence, *<sup>K</sup>p<sup>i</sup>* = **0**. Thus, *<sup>F</sup>*¯*Ki* contains only force/torque of virtual joints. Hence:

$$\begin{aligned} \bar{F}\_{\text{ext}} &= (J\_{\theta i}^{-T} \mathbf{K}\_{\theta i} I\_{\theta i}^{-1}) \Delta \bar{X}\_i \\ J\_{pi}^T \bar{F}\_{\text{ext}} &= \bar{0} \end{aligned} \tag{13}$$

The effect of passive joints is included in the stiffness model by inverting the following homogeneous relation matrix. This matrix is always invertible if det *JT θiJθ<sup>i</sup>* = 0.

$$
\begin{bmatrix}
(J\_{\theta i} \mathbf{K}\_{\theta i}^{-1} I\_{\theta i}^T) & J\_{pi} \\
J\_{pi}^T & \mathbf{0}
\end{bmatrix}^{-1} = \begin{bmatrix}
[\mathbf{K}\_{\zeta i}]\_{\theta \ge \theta} \vdots \underset{\theta \ge \theta}{\vdots} \cdots \\
\sim & \vdots \sim
\end{bmatrix} \tag{14}
$$

where *KCi* denotes stiffness matrix of *i*th in Cartesian space. Note that, *KCi* is rank deficient due to the passive joints. Cartesian stiffness matrix of the manipulator, *KC*, is computed as *KC* = ∑3 *i*=1 *KCi*. If external wrench is assumed to be *F*¯*ext* = 0, ¯ *KC* takes the following form.

$$
\mathbf{K}\_{\mathbb{C}} = \begin{bmatrix}
K\_{\mathbb{C}\_1}^{(11)} & 0 & 0 & 0 & K\_{\mathbb{C}\_1}^{(15)} & K\_{\mathbb{C}\_1}^{(16)} \\
0 & K\_{\mathbb{C}\_2}^{(22)} & 0 & K\_{\mathbb{C}\_2}^{(24)} & 0 & K\_{\mathbb{C}\_2}^{(26)} \\
0 & 0 & K\_{\mathbb{C}\_1}^{(33)} & K\_{\mathbb{C}\_3}^{(34)} & K\_{\mathbb{C}\_3}^{(35)} & 0 \\
0 & K\_{\mathbb{C}\_2}^{(24)} & K\_{\mathbb{C}\_2}^{(34)} & K\_{\mathbb{C}\_2}^{(44)} + K\_{\mathbb{C}\_3}^{(44)} & K\_{\mathbb{C}\_3}^{(45)} & K\_{\mathbb{C}\_2}^{(46)} \\
K\_{\mathbb{C}\_1}^{(15)} & 0 & K\_{\mathbb{C}\_3}^{(35)} & K\_{\mathbb{C}\_3}^{(35)} & K\_{\mathbb{C}\_1}^{(55)} + K\_{\mathbb{C}\_3}^{(56)} & K\_{\mathbb{C}\_1}^{(56)} \\
K\_{\mathbb{C}\_1}^{(16)} & K\_{\mathbb{C}\_2}^{(26)} & 0 & K\_{\mathbb{C}\_2}^{(46)} & K\_{\mathbb{C}\_1}^{(56)} & K\_{\mathbb{C}\_1}^{(66)} + K\_{\mathbb{C}\_2}^{(66)}
\end{bmatrix}
\tag{15}$$

For small deflections and loads, this matrix can be used without causing high errors. This matrix must be re-computed if |*F*¯*ext*| >> 0. *KC* is divided into 3 *×* 3 sub-matrices.

$$\mathbf{K}\_{\mathbf{C}} = \begin{bmatrix} \mathbf{K}\_{A} & \mathbf{K}\_{B} \\ \mathbf{K}\_{B}^{T} & \mathbf{K}\_{D} \end{bmatrix} \tag{16}$$

where *KA*, *KB*, *K<sup>T</sup> B*, and *KD* have the units of N/m, N/rad, N/rad, and Nm, respectively. The kinematic dimensions and material properties of the links are given in [47]. This study makes use of these parameters for stiffness evaluation of R-CUBE manipulator.

#### **3. Stiffness Performance Indices**

In this section, a review of the existing performance indices in literature is described, and a new performance index for isotropy evaluation of the stiffness matrix is proposed. Besides, usage of the indices for the R-CUBE manipulator is given in a table.

#### *3.1. Performance Indices in the Literature*

Stiffness performance indices are computed by obtaining the eigenvalues and eigenvectors. SVD operation may be used to obtain these properties as follows:

$$
\mathbb{K} = \mathbb{Q}\text{D}\mathbb{U} \tag{17}
$$

where *K* is an *n × n* stiffness matrix defined in Cartesian space for *n* = 1, 2, ..., 6. *Q* and *U* are orthogonal matrices whose columns are the eigenvectors of *K*. Since *K* is a symmetric matrix, *Q* = *UT*. *D* is a diagonal matrix whose elements denote the positive eigenvalues. These matrices are shown as:

$$\mathbf{Q} = \begin{bmatrix} \mathbb{E}\_1 & \mathbb{E}\_2 & \dots & \mathbb{E}\_6 \end{bmatrix} \tag{18}$$

$$\mathbf{D} = \text{diag}\left(\lambda\_1, \lambda\_2, \dots, \lambda\_6\right) \tag{19}$$

$$
\lambda\_1 \le \lambda\_2 \le \dots \le \lambda\_6 \tag{20}
$$

where *e*¯*i* and *λi* for *i* = 1, 2, ..., 6 denote the eigenvectors and eigenvalues, respectively. Regardless of the kinematic DoF of the manipulator, stiffness matrix is always a 6 *×* 6 matrix because compliant displacement may occur in any direction in 6 DoF space. Direct evaluation of a 6 *×* 6 stiffness matrix causes unit inconsistency in the results. A simple method is to use the normalized stiffness matrix [48]. The proposed method uses pre- and post-multiplication of the matrix with diagonal scaling matrices. Another method is introduced by Angeles [49]. He used the natural length to obtain normalized and unity matrix. Plücker Coordinates is preferred by Khan and Angeles to use dimensionally homogeneous space [50]. Thus, the dimensionally homogeneous matrix is obtained. A separate evaluation of rotation sensitive and translation sensitive matrices are introduced by Cardou et al. [51]. Hence, it is also possible to use interested sub-matrices of 6 *×* 6 stiffness matrix. In this way, a certain aspect of the stiffness model can be placed in focus. In our case, since we focus on the translational compliant displacements of the mobile platform, we used the top left corner of the stiffness matrix as indicated in Equation (16). Hence, in this sub-matrix, unit consistency is achieved. Nevertheless, this produces a lower dimensional stiffness matrix.

The graphical illustration of eigenvectors and eigenvalues generates an *n*-dimensional surface. This surface becomes a line for *n* = 1, an ellipse for *n* = 2, an ellipsoid for *n* = 3, and a hyper-ellipsoid for *n* = 4, 5, 6. The radii and axes of ellipsoids are defined by eigenvectors and eigenvalues as shown in Figure 3 for *n* = 3 [26].

**Figure 3.** Illustration of eigenvalues (*λi* for *i* = 1, 2, 3) and eigenvectors (*e*¯*i* for *i* = 1, 2, 3) as an ellipsoid for a 3 *×* 3 stiffness matrix.

Eigenvalues indicate the magnitude of the stiffness along their respective eigenvectors. Hence, the volume of this ellipsoid indicates the stiffness capacity. This volume is proportional to the determinant of the stiffness matrix [3,11,14,16,26–28,32]. Thus, stiffness value is commonly computed as:

$$\gamma\_V = \det(\mathbf{K}) = \prod\_{i=1}^n \lambda\_i \tag{21}$$

where *γV* denotes an average stiffness magnitude, and *λi* is the *i*th eigenvalue. *γV* equals to zero when inspected pose of the manipulator has at least one free motion direction. This case corresponds to a stiffness singularity. Thus, such a singular pose can be determined by investigating *γV*.

The Euclidean norm -.-*E*, which is also named as 2-norm -.-2, computes the square root of the largest positive eigenvalue (or singular value) of the square of a matrix, *KK<sup>T</sup>* shown as follows:

$$\|\mathbf{K}\|\_{E} = \max(\sqrt{\tilde{\lambda}\_i}) \tag{22}$$

where *λ* ˜ *i* denotes the *i*th eigenvalue of matrix *KKT*. This norm exhibits the maximum eigenvalue of *K* which is also the largest radius of the ellipsoid. The direction of this eigenvalue is the direction of the displacement which exhibits the highest stiffness. Euclidean norm of *K*−<sup>1</sup> reveals the minimum eigenvalue of *K*, which indicates the most compliant direction [14,17,27,29,30,33]. In optimization problems, it is common to only focus on minimum eigenvalue in the workspace in order to maximize it.

Other norms are 1 norm -.-1, infinity norm -.-∞ (also named as Chebyshev norm -.-*C*), and Frobenius norm -.-*F*. 1 norm computes the maximum of the summation of absolute values of column elements of a matrix. Infinity norm makes this computation for rows. Since the stiffness matrix is symmetric, both norms result in the same values. These norms denote a combined total resistive force and moment against a unit displacement in Cartesian space. Frobenius norm, on the other hand, has more meaning in terms of average stiffness. It focuses on diagonal elements of *KKT*. It is computed as follows:

$$\|\mathbf{K}\|\_{F} = \left(\sum\_{i=1}^{n} \sum\_{j=1}^{n} (\mathbf{K}\_{ij})^2\right)^{1/2} = \sqrt{\operatorname{tr}(\mathbf{K}\mathbf{K}^T)} = \sqrt{\lambda\_1^2 + \lambda\_2^2 + \dots + \lambda\_n^2} \tag{23}$$

where *tr* is the trace operator and subscript *ij* denotes the *ij*th element of the matrix. All norms result in higher values in stiff poses and lower values in compliant poses. Hence, they show similar relative distribution throughout the workspace. Therefore, an evaluation of one of the norms is sufficient to have an idea of stiffness distribution depending on manipulator pose.

Stiffness condition number, stiffness dexterity, or stiffness isotropy index is computed as the ratio of the maximum eigenvalue to the minimum eigenvalue [3,6,10,14,17,26,28,31]. This ratio reveals stiffness value distribution among eigenvectors. The minimum ratio of 1 indicates equal stiffness distribution. The calculation of this index is given as:

$$\gamma\_{\mathbb{C}} = \frac{\lambda\_n}{\lambda\_1} = \|\mathbb{K}\|\_E \left\|\mathbb{K}^{-1}\right\|\_E \tag{24}$$

where *γC* is stiffness condition number.

Another stiffness performance index is the uniformity index which compares the maximum and minimum values of *γC*. This index is formulated as follows:

$$\gamma\_{\rm CLI} = \frac{\max(\gamma\_{\rm C})}{\min(\gamma\_{\rm C})} \tag{25}$$

$$
\gamma\_{\rm Cu} \ge 1 \tag{26}
$$

where *γCU* is the uniformity index. Notice that, the minimum value of *γC* is bounded by 1 and so *γCU* ≥ 1. Here, *γCU* = 1 can only be achieved when the maximum and minimum values of *γC* are equal to each other. The only condition that satisfies this equality is when all *γC* in all poses are equal to each other.

The determination of both stiff and isotropic poses is a problem. A solution is proposed in [15] that evaluates the stiffness magnitude and isotropy, simultaneously. The formulation of this index is given below.

$$\gamma\_G = \frac{\lambda\_n^2 \lambda\_1^2}{\lambda\_n + \lambda\_1} \tag{27}$$

*γG* value increases as the manipulator approaches a stiff and isotropic pose. This index is more suitable for 2 *×* 2 stiffness matrices. For higher dimensions, this approach is not appropriate because the problem turns into a volumetric problem while *γG* solves a surface equation.

Energy index computation depends on whether a constant payload is applied or a constant compliant displacement is given to a manipulator. For a constant payload, stiff poses of the manipulator result in smaller deflections. Since the energy has a quadratic relationship with elastic displacements, stiff manipulators store less energy for the same payload. Therefore, the highest stored energy is observed when the elastic displacement is along the eigenvector direction that has the minimum eigenvalue for a payload. On the other hand, when a constant compliant displacement is given to manipulator, obviously stiff poses store more energy. It means the highest energy for a constant displacement is observed when this displacement is given along the eigenvector that has the highest eigenvalue. For a unit displacement, stored energy equals to the half of the maximum eigenvalue. The index formulation is given as:

$$\gamma\_c = \frac{1}{2} \Delta \mathbf{X}^T \mathbf{K} \Delta \mathbf{X} \quad \Rightarrow \gamma\_c = \frac{1}{2} \overline{e}\_n^T \mathbf{K} \ell\_n \quad \Rightarrow \gamma\_c = \frac{1}{2} \lambda\_n \quad \Rightarrow \gamma\_c = \frac{1}{2} ||\mathbf{K}||\_E \tag{28}$$

where *γ* is the energy index [5,7,11,52].

Unfortunately, all the above indices are computed numerically. In addition, they are local indices (except *γCU*). Hence, they are pose dependent, and they only indicate local performances. Thus, a globalization process is necessary for these indices [35]. This globalization can be achieved for the abovementioned indices to obtain an average value of performance of the whole workspace as follows:

$$
\bar{\mu} = \frac{\int \mu dv}{\int dv} \tag{29}
$$

$$
\sigma = \sqrt{\frac{\int (\mu - \bar{\mu})^2 dv}{\int dv}} \tag{30}
$$

where *μ* denotes any of the performance indices. *μ*¯ is an average value, and *σ* is the standard deviation. It is desired to have a low deviation to have a uniform stiffness performance distribution.

#### *3.2. Proposed Performance Index*

As can be seen in Section 3.1, many of the indices mostly rely on the maximum and minimum eigenvalues for performance evaluation. However, this causes a lack in proper stiffness isotropy evaluation for *n × n* dimensional stiffness matrix where *n* > 2. In this respect, a revision is required for stiffness isotropy index.

While stiffness condition number for a 2 *×* 2 stiffness matrix regarded as the ratios of radii of a stiffness ellipse, the mathematical meaning behind is the comparison of the areas of a circle and an ellipse. Area of this circle indicates an ideal (desired) performance value which has high isotropy and rigidity. Hence, the square of the maximum eigenvalue is an indication of the area of the circle. Area of the ellipse is the multiplication of all eigenvalues (determinant). The ratio between these areas gives the stiffness condition number. For a 2 *×* 2 stiffness matrix, this index is formulated as follows:

$$
\gamma\_{\mathbb{C}} = \frac{\lambda\_2^2}{\lambda\_1 \lambda\_2} = \frac{\lambda\_2}{\lambda\_1} \tag{31}
$$

$$
\gamma c \ge 1 \tag{32}
$$

The current *γC* is a 2 DoF evaluation approach and has a lack of performance in the evaluation of 3 or more DoF problems. Therefore, *γC* cannot distinguish the performances of poses whose eigenvalues *λ*1 = *λ*2 < *λ*3, *λ*1 < *λ*2 = *λ*3, or *λ*1 < *λ*2 < *λ*3 for 3 *×* 3 matrix. Poses that have closer eigenvalues to maximum one are more isotropic. Hence, *γC* must be revised for higher DoF evaluation. Higher dimensional isotropy index is the comparison of volumes of an ideal *n*-dimensional sphere and an *n*-dimensional ellipsoid. In this regard, the extension of the condition number for *n × n* stiffness matrix is a volumetric condition number or volumetric isotropy index that is formulated as follows:

$$\gamma\_I = \frac{\|\mathbf{K}\|\_E^n}{\det(\mathbf{K})} \tag{33}$$

$$
\gamma\_1 \ge 1 \tag{34}
$$

where *γI* volumetric stiffness isotropy index and *n* is the dimension of the stiffness matrix, *n*th power of -*K*-*E* indicates the volume of the ideal sphere while det(*K*) represented the volume of the ellipsoid. *γI* = 1 indicates the isotropic poses in terms of stiffness. Since this index also considers the intermediate eigenvalue contribution, the ambiguity of stiffness performance between the poses that have same *γC* can be accurately distinguished. Hence, it is expected to have a lower standard deviation of *σI* compared to *σC*. The standard deviation in such a comparison can be regarded as the preciseness of the performance index.

A volumetric uniformity index *γIU* may also be defined as an extension to *γI*. The uniformity index is defined as follows:

$$\gamma\_{III} = \frac{\max(\gamma\_I)}{\min(\gamma\_I)}\tag{35}$$

$$\gamma\_{\rm mf} \ge 1\tag{36}$$

Since the preciseness of the volumetric stiffness isotropy index is expected to be relatively better, the volumetric uniformity index is expected to be more accurate compared to *γCU*.

#### *3.3. Construction of Performance Indices for R-CUBE*

Previously mentioned indices are utilized for the stiffness evaluation R-CUBE manipulator, and they are tabulated in Table 1. Since the mobile platform of the manipulator has only translational DoF, translational compliant displacements are our primary interest. Hence, *KA* sub-matrix is in the focal point of this study.


**Table 1.** Utilized stiffness performance indices.

#### **4. Results and Discussion**

In this section, only the results on isotropy index calculations are given to present a comparison between the conventional and proposed isotropy index. *γI*, and *γC* values are computed throughout the workspace, and their results are compared with each other. These computations are carried out by using *KA* sub-matrix. A normal distribution plot is presented for both indices to illustrate the relative preciseness of the proposed index *γI*. Isotropy indices are computed for each discrete pose of the workspace. Then, the computed indices are illustrated via color mapping.

The maximum and minimum values, average values, standard deviations of isotropy indices are presented in Table 2. Having the minimum value of 1 for *γC* and *γI* indices shows that there at least one isotropic pose. By definition of both indices, the isotropy is observed in the same poses. Even though both indices compute the same stiffness matrix, their maximum values differ from each other as expected. The effect of intermediate eigenvalue causes this difference. These maximum values indicate the worst isotropy performance. The number of worst isotropic poses with respect to *γC* is higher compared to *γI* because *γC* does not consider the intermediate eigenvalues. However, the worst isotropic poses concerning *γI* are also the worst isotropic poses with respect to *γC*, but vice-versa is not necessarily valid for all least isotropic poses.

The values of *γC* and *γI* are normalized as *γ*∗ *C* and *γ*∗ *I* to carry out a fair comparison among them. Accordingly, *μ*¯∗ *I* , *μ*¯∗ *C*, *σ*<sup>∗</sup> *I* , and *σ*<sup>∗</sup> *C* are obtained in a normalized space. *γC* and *γI* are normalized such that the worst performance is denoted by 1 and the highest isotropy is given by 0.


**Table 2.** Results of isotropy indices.

The normal distribution of both indices is compared in a normalized space and shown in Figure 4. Probability density in the vertical axis indicates the number of poses that have the isotropy value denoted in the horizontal axis. Mean values *μ*∗ *I* , and *μ*∗ *C* indicates the average isotropy of the workspace. *γ*∗ *I* results indicate that the workspace is more isotropic with respect to *γ*∗ *C* results since *μ*∗ *I* < *μ*∗ *C*. In addition, the probability distribution of isotropic poses is higher by using *γI* index. The standard deviation *σ*<sup>∗</sup> *I* is lower than *σ*<sup>∗</sup> *C*. Hence, *γ*∗ *I* is more sensitive to isotropy changes between poses, and this new index can detect slightest changes.

**Figure 4.** Normal distribution of *γI* and *γC*.

For illustration purposes, 9 planes are selected in the workspace. These planes are parallel to*u*1 −*u*2, *u*2 −*u*3, and*u*1 −*u*3 orthogonal planes. They are divided into 3 groups. In each group, 3 planes are intersect at [−60,−60,−<sup>60</sup>] mm, [0, 0, 0] mm, [60, 60, 60] mm locations of the workspace. Figures 5a and 6a, show the outermost surfaces of the workspace where the planes are placed at [60, 60, 60] mm. Figures 5b and 6b, show the planes that are coincident at [0, 0, 0] mm. Figures 5c and 6c, show the most inner surfaces of the workapce at [−60,−60,−<sup>60</sup>] mm. Due to the symmetric topology of the manipulator, computed performance indices have the same value in several poses. Hence, a symmetric distribution is observed in the results of both indices throughout the workspace.

Due to the definition of both indices, lower values (represented with blue zones) indicate higher isotropy. These areas mainly are observed at the innermost corner, the outermost corner, and at the center of the workspace. However, isotropy distribution differs for *γ*∗*I* and *<sup>γ</sup>*<sup>∗</sup>*C* for the rest of the workspace. The transition between low isotropy and high isotropy is more smooth in *<sup>γ</sup>*<sup>∗</sup>*C* compared to *γ*∗*I* . This smooth transition is a consequence of relatively less precise results obtained with *γC*. Hence, *γC* cannot precisely distinguish the isotropy levels among poses. However, the transition from low to high isotropy in *γ*∗*I*results illustration is sharper and more distinguishable compared to *<sup>γ</sup>*<sup>∗</sup>*C*.

**Figure 5.** Isotropy distribution: (**a**) Outer surfaces, (**b**) middle surfaces, (**c**) inner surfaces.

**Figure 6.** Volumetric isotropy distribution: (**a**) Outer surfaces, (**b**) middle surfaces, (**c**) inner surfaces.

Figures 7 and 8 show the iso-curves in *u*1 −*u*2 planes for *<sup>γ</sup>*<sup>∗</sup>*C*, and *γ*∗*I* , respectively. These planes are positioned along *u*3 and they are located at 60 mm,0 mm,−60 mm. Figure 8 shows that the covered area by the iso-curves for *γ*∗*I* = 0.1, 0.2 is higher than the area of *<sup>γ</sup>*<sup>∗</sup>*C* = 0.1, 0.2 in Figure 7. Hence, *γI* captures more isotropic poses. Accordingly, iso-curves for *<sup>γ</sup>*<sup>∗</sup>*C* = 0.3, 0.4, ..., 0.8 cover more area than the area covered by the same value of *γ*∗*I* . This indicates that the manipulator is less isotropic with respect to *γC*. In addition, notice that the distance between the iso-curves are less for *γ*∗*I* in Figure 8 compared to *<sup>γ</sup>*<sup>∗</sup>*C*in Figure 7. Therefore, *γI* is more sensitive to isotropy changes.

**Figure 7.** Normalized isotropy iso-curves: (**a**) 60 mm along *u*3, (**b**) 0 mm along *u*3, (**c**) −60 mm along *u*3.

**Figure 8.** Normalized volumetric isotropy iso-curves: (**a**) 60 mm along *u*3, (**b**) 0 mm along *u*3, (**c**) −60 mm along *u*3.
