2. Structural design

As mentioned in previous study of the authors [Reference Gezgin, Özbek, Güzin, Ağbaş and Gezer19], kinematic structure of the proposed manipulator includes dual concentric arcs as the input links that are connected to the fixed ground orthogonally by utilizing revolute pairs. Each input link houses a single platform that can slide along the circumference of an arc and both of the platforms are assembled together by another revolute joint whose axis passes through the isocenter of the system. At this isocenter, all of the axes of the joints of the manipulator are intersected and the center of the arcs resides (Fig. 1).

Figure 1. Kinematic representation of proposed spherical manipulator.

Proposed system mainly designed to provide proper and precise positioning of the needle guide during brain biopsy operations by considering tumor locations inside the brain. As the manipulator has a spherical structure, once its isocenter gets aligned with the central volume that resides at the intersection of sagittal, coronal and transverse planes of the human skull, required alignment of the needle guide with selected tumor can easily be acquired by using inverse task equations of the system.

Although the first version of the system was proved to be capable of carrying out positioning tasks in a brain biopsy mock-up trials [Reference Gezgin, Özbek, Güzin, Ağbaş and Gezer19], its conceptual design has been updated to increase versatility of the system. Thus, in addition to brain biopsy, it can now be utilized on various brain surgery operations such as craniotomy, neuroendoscopy, and deep brain stimulation.

As seen in Fig. 2, not only link designs of the manipulator was improved to increase reachable workspace and manufacturability of the system but also its platform was conceptually redesigned to allow attachment of different tools in order to adapt various brain surgery applications by utilizing additional degree of freedom along the normal axis passing through the isocenter. Thus workspace of the system has also been improved from a spherical surface to a spherical volume (Fig. 3).

Figure 2. Concept of the updated design and platform module that can be integrated to the upper platform with different tools.

Figure 3. Workspace of the manipulator with and without the attached module.

Figure 4. Percentages of possible tumor locations and structural design.

Structural dimensioning of the new design was carried out by considering percentages of tumor locations inside human brain.

Gould [Reference Gould20] showed that malignant primary brain tumor locations are mostly concentrated on the frontal lobe (26%) of the brain. Thus during this study frontal lobe was tried to be covered as much as possible by avoiding a bulk design with large dimensions. Frontal lobe coverage also allows for the manipulator to reach parietal lobe (12%), ventricle (2%), and some part of the temporal lobe (19%). In order to contain adult human head inside the manipulator workspace by considering important tumor locations, lower and upper arc radiuses of the manipulator should be properly decided. In light of this, anthropometric survey of US army personnel [Reference Gordon, Churchill, Clauser, Bradtmiller, McConville, Tebbetts and Walker21] was utilized, where more than 1700 men and 2000 women sampled for measurements. As represented in Table I, averages of given measurements were taken into consideration and radiuses of the lower and upper circular arcs were chosen as 125 and 160 mm, respectively. Since the averages of male anthropometric measurements are bigger than the averages of female anthropometric measurements, average male data sets were used in the design.

Table I. Anthropometric data of US army personel (mm) [Reference Gordon, Churchill, Clauser, Bradtmiller, McConville, Tebbetts and Walker21].

In terms of utilization during surgical procedures, the vital difference between earlier prototype and its proposed enhanced version is the fact that the latter one has the potential to carry out regular surgical procedures such as opening a burr hole, creating a bone flap, inserting an electrode or a biopsy needle to the brain while the former one just assist the surgeon to have a proper and precise alignment with respect to the target. From this point of view, while it might be enough for the former version to be controlled by considering only position and orientation of the platform, enhanced modular version should be controlled by taking additional considerations in terms of forces and torques in order to increase safety and reliability. Thus, it is very crucial to perform dynamic analysis of the proposed manipulator system. In light of this, dynamic analysis of the parallel section of the manipulator was carried out by using both simulation environment and theoretical approach during the designed case study of opening a bone flap on human skull.

3. Modeling in simulation environment

Although addition of the platform modules fundamentally shifts manipulators parallel structure to a hybrid (parallel and serial) one, they were assumed as platform payloads for the selected case study. Thus, dynamic analysis of the system was realized by only focusing on the 2-DOF parallel section of the manipulator that is carrying a payload with its top platform (Fig. 5).

Figure 5. Simplification of the structure for the dynamic analysis.

Prior to the theoretical approach, for the ease of use, dynamic model of the system was first carried out in a MATLAB environment. Once the updated 3D CAD model of the system was prepared with proper assembly constraints, it was imported into SimMechanics module of the software (Fig. 6) and trajectory planning of the case study was performed.

Figure 6. Modelled manipulator at MATLAB SimMechanics environment.

Figure 7. Landmark points and their projections on upper platform workspace.

3.1. Case study and trajectory planning

For the aim of extracting torque data from the actuators of the system in case of a given platform trajectory, it was decided to construct a case study that involves surgical application of opening a bone flap on human skull. As seen in Fig. 7, four desired landmark points were selected to generate a simple bone flap on a human skull model within the workspace of the manipulator. In order to have the end effector of the system to cut desired bone section carefully with respect to the landmark constraints, trajectory of the upper platform should be created with respect to the given constraints. Thus prior to the trajectory generation, selected landmark points (A₁, A₂, A₃, and A₄) on the skull were projected to the spherical workspace of the upper platform by using Eq. (1) as

(1) \begin{equation} \mathrm{B}_{ix}=\mathrm{A}_{ix}\frac{\left| \mathbf{r}_{u}\right| }{\left| \mathbf{r}_{\mathrm{Ai}}\right| },\,\,\mathrm{B}_{iy}=\mathrm{A}_{iy}\frac{\left| \mathbf{r}_{u}\right| }{\left| \mathbf{r}_{\mathrm{Ai}}\right| },\,\,\,\mathrm{B}_{iz}=\mathrm{A}_{iz}\frac{\left| \mathbf{r}_{u}\right| }{\left| \mathbf{r}_{\mathrm{Ai}}\right| }\quad i\rightarrow 1,2,3,4 \end{equation}

where A _ix , A _iy , A _iz , B _ix , B _iy , and B _iz are the x, y, z coordinates of i ^th A and B points, r_Ai is the distance between i ^th A point and isocenter of the system and r _u is the radius of spherical workspace of the upper platform that will be measured up to the center of mass of the platform from the isocenter throughout the study due to the symmetry.

Table II. Boundary conditions for the given trajectory.

Due to the fact that upper platform should pass from points B₁, B₂, B_3, and B₄, in order to generate necessary boundary conditions, required actuator positions to reach these points should also be calculated by using kinematic representation (Fig. 1) and inverse task equations of the manipulator that were explained in detail in ref. [Reference Gezgin, Özbek, Güzin, Ağbaş and Gezer19] (Eq. (2)).

(2) \begin{equation} \theta _{1i}=\tan ^{-1}\frac{-\mathrm{B}_{iy}}{\mathrm{B}_{iz}},\,\,\,\,\,\theta _{2i}=\tan ^{-1}\frac{\mathrm{B}_{ix}}{\mathrm{B}_{iz}}\quad i\rightarrow 1,2,3,4 \end{equation}

During the cutting process, upper platform of the manipulator will start its motion from the initial position B₁ at time t ₁, pass through via points B₂, B₃, B₄ and it will return back to its initial position by coming to halt at time t ₅. Thus, the trajectory of the manipulator will include four separate segments as B₁-B₂, B₂-B₃, B₃-B₄, and B₄-B₁. In light of this, to find proper joint angle functions in joint space, boundary conditions shown in Table II were selected and applied to solve the unknowns of cubic polynomial functions as

(3) \begin{equation} \theta _{1j}(t)=a_{j}t^{3}+b_{j}t^{2}+c_{j}t+d_{j},\,\,\,\,\theta _{2j}(t)=e_{j}t^{3}+f_{j}t^{2}+g_{j}t+h_{j}\,\,\,\,\,j\rightarrow 1,2,3,4 \end{equation}

where j represents the segment number and each function is the actuator joint angle function on j ^th segment.

It is clear that Eq. (3) has 32 unknowns. Thus, utilizing given boundary conditions, these unknowns can be solved (Eq. (4)) by forming 32 polynomial equations to reach joint angle functions for each individual segment. It should be noted that, at each intermediate via points, not only positions but also velocities and accelerations are equalized in order to have continuity throughout the functions.

(4) \begin{equation} \underset{32x32}{\mathbf{D}}\underset{32x1}{\left[\begin{array}{l} a_{1}\\ .\\ .\\ h_{4} \end{array}\right]}=\underset{32x1}{\left[\begin{array}{l} \theta _{11}\\ .\\ .\\ \theta _{24} \end{array}\right]},\,\,\,\,\,\underset{32x1}{\left[\begin{array}{l} a_{1}\\ .\\ .\\ h_{4} \end{array}\right]}=\underset{32x32}{\mathbf{D}^{-1}}\underset{32x1}{\left[\begin{array}{l} \theta _{11}\\ .\\ .\\ \theta _{24} \end{array}\right]} \end{equation}

For the sake of acquiring numerical values, Eq. (4) was solved by utilizing parameters of the created bone flap geometry example considering 10 s time intervals between each segment (Table III).

Table III. Simulation constraints and solved parameters.

Using calculated parameters in conjunction with Eq. (3), joint angle functions and joint rates were drawn in separate graphs for each actuator (Fig. 8a). Moreover, by the help of systems direct task equations (Eq. (5)), platform trajectory and its projection on skull model were also demonstrated (Fig. 8b).

(5) \begin{equation} \mathrm{B}_{x}=\frac{r_{u}C_{1}S_{2}}{\sqrt{{C_{2}}^{2}+{C_{1}}^{2}{S_{2}}^{2}}},\,\,\mathrm{B}_{y}=-\frac{r_{u}S_{1}C_{2}}{\sqrt{{C_{2}}^{2}+{C_{1}}^{2}{S_{2}}^{2}}},\,\,\mathrm{B}_{z}=\frac{r_{u}C_{1}C_{2}}{\sqrt{{C_{2}}^{2}+{C_{1}}^{2}{S_{2}}^{2}}},\,\,C_{i}=\cos \theta _{i},\,\,S_{i}=\sin \theta _{i} \end{equation}

4. Dynamic analysis of the manipulator by using Lagrange method

In order to verify acquired simulation results for the given case study, theoretical dynamic analysis of the manipulator by utilizing Lagrangian formulation was carried out. Lagrangian formulation of the first kind that is written by a set of redundant coordinates can be seen in Eq. (6),

(6) \begin{equation} \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_{j}}\right)-\frac{\partial L}{\partial q_{j}}=Q_{j}+{\sum }_{i=1}^{k}\lambda _{i}\frac{\partial \unicode[Times]{x03B3} _{i}}{\partial q_{j}},\,\,\,\,\,L=K-U,\,\,\,\,\,i\rightarrow 1,\ldots,k,\,\,\,\,\,j\rightarrow 1,\ldots,n \end{equation}

where q _j is the j ^th generalized coordinate, Q _j is the actuator force or torque, λ _i is the i ^th Lagrangian multiplier, Γ_i is the i ^th constraint function, K and U are the kinetic and potential energies of the manipulator, respectively, k is the number constraint functions and n is the total number of generalized coordinates. As mentioned by Tsai [Reference Tsai4], for the ease of use, Eq. (6) can be arranged into two sets of equations including k equations related with redundant coordinates and n-k equations related with actuated joint variables, as shown in Eq. (7),

(7) \begin{equation} \begin{array}{c} {\sum }_{i=1}^{k}\lambda _{i}\dfrac{\partial \unicode[Times]{x03B3} _{i}}{\partial q_{j}}=\dfrac{d}{dt}\left(\dfrac{\partial L}{\partial \dot{q}_{j}}\right)-\dfrac{\partial L}{\partial q_{j}}-\hat{Q}_{j},\,\,\,\,\,\,i\rightarrow 1,\ldots,k,\,\,\,\,\,j\rightarrow 1,\ldots,k\\[14pt] Q_{j}=\dfrac{d}{dt}\left(\dfrac{\partial L}{\partial \dot{q}_{j}}\right)-\dfrac{\partial L}{\partial q_{j}}-{\sum }_{i=1}^{k}\lambda _{i}\dfrac{\partial \unicode[Times]{x03B3} _{i}}{\partial q_{j}},\,\,\,\,\,\,i\rightarrow 1,\ldots,k,\,\,\,\,\,j\rightarrow k+1,\ldots,n \end{array} \end{equation}

where $\hat{Q}_{j}$ represents generalized force component due to an external force, and it will be assumed to be zero for this preliminary study. As seen in Eq. (7), for the aim of determining necessary actuator forces or torques, constraint equations of the manipulator system should be derived by the help of manipulator kinematics. Also Lagrangian multipliers should be found as a follow-up.

As seen in Fig. 10, proposed manipulator system is based on a spherical five bar linkage, where all of the twist angles are 90°. Prior to the derivation of constraint equations, generalized coordinates of the manipulator system were selected as three redundant coordinates of the upper platform position B _x , B _y , and B _z along with two actuated joint parameters θ ₁ and θ ₂. Due to the fact that three redundant coordinates were selected to simplify Lagrangian dynamics, the same number of constraint equations should be derived.

It is clear that upper platform of the manipulator always move on a spherical surface with a radius r due to the manipulators spherical structure. As a result normal distance between the isocenter of the system and the platform mass center will always becomes

(8) \begin{equation} \sqrt{{\mathrm{B}_{x}}^{2}+{\mathrm{B}_{y}}^{2}+{\mathrm{B}_{z}}^{2}}=r_{u} \end{equation}

Figure 8. (a) Trajectory of the upper platform and its projection on human skull model, (b) joint angle functions and joint rates.

Figure 9. (a) SimMechanics model with joint angle functions as inputs, (b) necessary actuator torque values that are required to manipulate upper platform.

Figure 10. (a) Spherical five bar linkage, (b) derivation of constraint equations.

Thus, using Eq. (8), the first constraint equation of the system can be written as

(9) \begin{equation} \unicode[Times]{x03B3} _{1}\rightarrow {\mathrm{B}_{x}}^{2}+{\mathrm{B}_{y}}^{2}+{\mathrm{B}_{z}}^{2}-{r_{u}}^{2}=0 \end{equation}

Before proceeding forward from this point, input links should clearly be focused on. Note that, although they are assembled together by a revolute joint, upper and lower platform of the manipulator slides on the surfaces of separate semi-circular links. During these motions, input link on x axis related with the upper platform is actuated by an angle θ ₁, while the input link on y axis related with the lover platform is actuated by an angle θ ₂.

In order to find the second constraint equation, let us consider that the arc of the upper platform is actuated first from the manipulators initial configuration (both of the input links stay perpendicular to the xy plane). If the link rotates around the x axis by θ ₁ , the normal vector of the arc plane N _u will also carry out the same rotation. This action actually slides the lower platform of the system on its own input link. If the second inputs actuation by θ ₂ occurs after this motion, the upper platform will reach its final position by the rotation around the normal N _u . Thus, as seen in Fig. 10, normal vector N _u will always constrained to be perpendicular to the position of the upper platform. In its initial configuration, unit vector of the normal N _u lies on the positive y axis as

(10) \begin{equation} \hat{\mathbf{N}}_{\mathbf{u}}=\mathbf{j} \end{equation}

If it is rotated around the x axis along with the input link of the upper platform by the angle θ ₁, its final position can be calculated by using a quaternion operator q( )q as

(11) \begin{equation} \hat{\mathbf{N}}_{\mathbf{u}\mathbf{2}}=\mathbf{q}\left(\hat{\mathbf{N}}_{\mathbf{u}}\right)\mathbf{q}^{\mathbf{-}\mathbf{1}}=\cos \theta _{1}\mathbf{j}+\sin \theta _{1}\mathbf{k} \end{equation}

where, $\mathbf{q}=\cos \dfrac{\theta _{1}}{2}+\mathbf{i}\sin \dfrac{\theta _{1}}{2}$ , and $\mathbf{q}^{\mathbf{-}\mathbf{1}}=\cos \dfrac{\theta _{1}}{2}-\mathbf{i}\sin \dfrac{\theta _{1}}{2}$ . If the result of Eq. (11) and unit vector of the platforms position are multiplied by dot product, due to their orthogonality the result will always be zero.

(12) \begin{equation} \hat{\mathbf{N}}_{\mathbf{u}\mathbf{2}}.\hat{\mathbf{r}}=\left(\cos \theta _{1}\mathbf{j}+\sin \theta _{1}\mathbf{k}\right).\frac{(\mathrm{B}_{x}\mathbf{i}+\mathrm{B}_{y}\mathbf{j}+\mathrm{B}_{z}\mathbf{k})}{r_{u}}=\frac{\mathrm{B}_{y}\cos \theta _{1}+\mathrm{B}_{z}\sin \theta _{1}}{r_{u}}=0 \end{equation}

Figure 11. Component separation.

Using Eq. (12), the second constraint equation of the system can be written as

(13) \begin{equation} \unicode[Times]{x03B3} _{2}\rightarrow \mathrm{B}_{y}\cos \theta _{1}+\mathrm{B}_{z}\sin \theta _{1}=0 \end{equation}

Similarly if the same procedure will be applied in reverse order, such as rotating the second input link by θ ₂ followed by the rotation of the other input link by θ ₁, final constraint equation of the system can also be derived as,

(14) \begin{equation} \begin{array}{c} \hat{\mathbf{N}}_{\mathbf{L}}=\mathbf{i}\\ \hat{\mathbf{N}}_{\mathbf{L}\mathbf{2}}=\mathbf{q}\left(\hat{\mathbf{N}}_{\mathbf{L}}\right)\mathbf{q}^{\mathbf{-}\mathbf{1}}=\cos \theta _{2}\mathbf{i}-\sin \theta _{2}\mathbf{k}\\ \unicode[Times]{x03B3} _{3}\rightarrow \mathrm{B}_{x}\cos \theta _{2}-\mathrm{B}_{z}\sin \theta _{2}=0 \end{array} \end{equation}

where $\hat{\mathbf{N}}_{\mathbf{L}}$ is the unit vector of the arc plane normal related with the lover platform and in its initial configuration, it lies on the positive x axis (Fig. 10).

Figure 12. Final floating coordinate frames of the upper and the lower platforms.

Figure 13. (a) Necessary actuator torque values by Lagrange method that are required to manipulate upper platform, (b) comparison between SimMechanics model and Lagrange method as difference graphs.

Prior to the utilization of Eq. (7), manipulator system was separated into four main components as upper platform, lower platform, and and two input links with their shafts assembled for the ease of use (Fig. 11).

In order to form Lagrangian formulation, kinetic energies of the manipulator for each of the separated components were written as

\begin{equation*} \begin{array}{c} K=K_{u}+K_{L}+K_{ua}+K_{La}\\[12pt] K_{u}=\dfrac{1}{2}m_{u}\left({\dot{\mathrm{B}}_{x}}^{2}+{\dot{\mathrm{B}}_{y}}^{2}+{\dot{\mathrm{B}}_{z}}^{2}\right)+\dfrac{1}{2}{\mathbf{w}_{u}}^{\mathrm{T}}\left({}_{u}^{0}{\mathbf{R}}{}\mathbf{I}_{u}\left({}_{u}^{0}{\mathbf{R}}{}\right)^{\mathrm{T}}\right)\mathbf{w}_{u}\\[12pt] K_{L}=\dfrac{1}{2}m_{L}\left({\dot{\mathrm{D}}_{x}}^{2}+{\dot{\mathrm{D}}_{y}}^{2}+{\dot{\mathrm{D}}_{z}}^{2}\right)+\dfrac{1}{2}{\mathbf{w}_{L}}^{\mathrm{T}}\left({}_{L}^{0}{\mathbf{R}}{}\mathbf{I}_{L}\left({}_{L}^{0}{\mathbf{R}}{}\right)^{\mathrm{T}}\right)\mathbf{w}_{L} \end{array} \end{equation*}

(15) \begin{equation} \begin{array}{c} K_{ua}=\dfrac{1}{2}m_{ua}\left(\left| \mathbf{w}_{ua}\right| r_{ua}\right)^{2}+\dfrac{1}{2}{\mathbf{w}_{ua}}^{\mathrm{T}}\left({}_{ua}^{0}{\mathbf{R}}{}\mathbf{I}_{ua}\left({}_{ua}^{0}{\mathbf{R}}{}\right)^{\mathrm{T}}\right)\mathbf{w}_{ua}\\[12pt] K_{La}=\dfrac{1}{2}m_{La}\left(\left| \mathbf{w}_{La}\right| r_{La}\right)^{2}+\dfrac{1}{2}{\mathbf{w}_{La}}^{\mathrm{T}}\left({}_{La}^{0}{\mathbf{R}}{}\mathbf{I}_{La}\left({}_{La}^{0}{\mathbf{R}}{}\right)^{\mathrm{T}}\right)\mathbf{w}_{La} \end{array} \end{equation}

where K _u , K _L , K _ua , and K _La are the kinetic energies, m _u , m _L , m _ua , and m _La are the masses, w _u , w _L , w _ua , and w _La are the angular velocity vectors, I _u , I _L , I _ua , and I _La are tensor of inertias about the center of masses, and ${}_{u}^{0}{\mathbf{R}}{},\,\,{}_{L}^{0}{\mathbf{R}}{},\,\,{}_{ua}^{0}{\mathbf{R}}{},\mathrm{and}\,\,{}_{La}^{0}{\mathbf{R}}{}$ are the rotation matrices with respect to the base frame of upper platform, lower platform, upper input link, and lower input link, respectively. Also, D_x, D_y, and D_z are the coordinates of the lower platforms mass center r _ua and r _La are the radiuses of the upper and lower input links mass centers, respectively, measured from the isocenter.

Among the parameters of Eq. (15), $\mathbf{w}_{ua},\,\mathbf{w}_{La},\,{}_{ua}^{0}{\mathbf{R}}{},\,\mathrm{and}\,{}_{La}^{0}{\mathbf{R}}{}$ are the easiest ones to be represented as they are related with sole rotations around x axis for the upper arc and y axis for the lower arc. Thus, these parameters can easily be represented in matrix form as

(16) \begin{equation} \mathbf{w}_{ua}=\left[\begin{array}{l} \dot{\theta }_{1}\\ 0\\ 0 \end{array}\right],\,\,\mathbf{w}_{La}=\left[\begin{array}{l} 0\\ \dot{\theta }_{2}\\ 0 \end{array}\right],\,\,{}_{ua}^{0}{\mathbf{R}}{}=\left[\begin{array}{l@{\quad}l@{\quad}l} 1 & 0 & 0\\ 0 & C_{1} & -S_{1}\\ 0 & S_{1} & C_{1} \end{array}\right],\,\,\,{}_{La}^{0}{\mathbf{R}}{}=\left[\begin{array}{l@{\quad}l@{\quad}l} C_{2} & 0 & S_{2}\\ 0 & 1 & 0\\ -S_{2} & 0 & C_{2} \end{array}\right] \end{equation}

On the other hand $\mathbf{w}_{u},\,\mathbf{w}_{L},\,\,{}_{u}^{0}{\mathbf{R}}{},\mathrm{and}\,\,{}_{L}^{0}{\mathbf{R}}{}$ should be represented in a more complicated way. In Fig. 12, floating coordinate frames of upper and lower platforms at their final positons are represented following the input rotations of θ ₁ and θ ₂.

Table IV. Deviation between torques and differences.

Figure 14. Change of friction forces during given trajectory.

Figure 15. Trajectory of the joints during hardware verification.

Figure 16. First manufactured prototype of the spherical manipulator.

Figure 17. Teflon guides on rails and needle bearing connecting top and bottom platforms.

In order to represent ${}_{u}^{0}{\mathbf{R}}{},\mathrm{and}\,\,{}_{L}^{0}{\mathbf{R}}{}$ , unit vectors of the floating coordinate frames should be projected on to the base frame. Due to the fact that representations of $\hat{\mathbf{z}}_{u}$ and $\hat{\mathbf{z}}_{L}$ unit vectors of the floating frames in base frame are always the unit vector of platform position vector, their projections can be instantly written as its x, y, and z components. For the upper platform $\hat{\mathbf{y}}_{u}$ unit vector lie on the yz plane by an angle θ ₁ measured from the y axis and for the lower platform $\hat{\mathbf{x}}_{L}$ unit vector lie on the xz plane by an angle θ ₂ measured from the x axis. Using this information it is also easy to project them to the base frame coordinates. Lastly, remaining projections of $\hat{\mathbf{x}}_{u}$ and $\hat{\mathbf{y}}_{L}$ can be calculated by the cross products of projected $\hat{\mathbf{y}}_{u},\,\,\hat{\mathbf{z}}_{u}$ , and $\hat{\mathbf{z}}_{L},\,\,\hat{\mathbf{x}}_{L}$ , respectively, to find final rotation matrices as

(17) \begin{equation} \,{}_{u}^{0}{\mathbf{R}}{}=\left[\begin{array}{l@{\quad}l@{\quad}l} \dfrac{C_{1}B_{z}-S_{1}B_{y}}{r_{u}} & 0 & \dfrac{B_{x}}{r_{u}}\\[13pt] \dfrac{S_{1}B_{x}}{r_{u}} & C_{1} & \dfrac{B_{y}}{r_{u}}\\[13pt] \dfrac{-C_{1}B_{x}}{r_{u}} & S_{1} & \dfrac{B_{z}}{r_{u}} \end{array}\right],\,\,\,{}_{L}^{0}{\mathbf{R}}{}=\left[\begin{array}{l@{\quad}l@{\quad}l} C_{2} & \dfrac{-S_{2}B_{y}}{r_{u}} & \dfrac{B_{x}}{r_{u}}\\[13pt] 0 & \dfrac{C_{2}B_{z}+S_{2}B_{x}}{r_{u}} & \dfrac{B_{y}}{r_{u}}\\[13pt] -S_{2} & \dfrac{-C_{2}B_{y}}{r_{u}} & \dfrac{B_{z}}{r_{u}} \end{array}\right] \end{equation}

If the derivatives of the rotation matrices shown in Eq. (17) are multiplied by their inverses, skew symmetric matrices of angular velocities can also be found.

(18) \begin{equation} \begin{array}{c} \,{}_{u}^{0}{\dot{\mathbf{R}}}{}{}_{u}^{0}{\mathbf{R}}{_{}^{-1}}=\left[\begin{array}{l@{\quad}l@{\quad}l} 0 & -w_{uz} & w_{uy}\\[5pt] 0 & 0 & -w_{ux}\\[5pt] -w_{uy} & w_{ux} & 0 \end{array}\right],\,\,\,{}_{L}^{0}{\dot{\mathbf{R}}}{}{}_{L}^{0}{\mathbf{R}}{_{}^{-1}}=\left[\begin{array}{l@{\quad}l@{\quad}l} 0 & -w_{Lz} & w_{Ly}\\[5pt] 0 & 0 & -w_{Lx}\\[5pt] -w_{Ly} & w_{Lx} & 0 \end{array}\right]\\[10pt] \mathbf{w}_{u}=\left[\begin{array}{l} w_{ux}\\[5pt] w_{uy}\\[5pt] w_{uz} \end{array}\right],\,\,\,\mathbf{w}_{L}=\left[\begin{array}{l} w_{Lx}\\[5pt] w_{Ly}\\[5pt] w_{Lz} \end{array}\right] \end{array} \end{equation}

After the parameters of Eq. (15) were derived, potential energies of the separated components can also be written as

(19) \begin{equation} \begin{array}{c} U=U_{u}+U_{L}+U_{ua}+U_{La}\\[5pt] U_{u}=m_{u}g\,\mathrm{B}_{z},\,\,\,U_{L}=m_{L}g\,D_{z},\,\,\,U_{ua}=m_{ua}g\,r_{ua}{C}_{1},\,\,\,U_{La}=m_{La}g\,r_{La}{C}_{2} \end{array} \end{equation}

where g is the gravitational acceleration along negative z axis and U _u, U _L , U _ua , and U _La are the potential energies of upper platform, lower platform, upper input link, and lower input link, respectively. Substituting Eqs. (15) and (19) to the Eq. (7), actuator torques can be derived after calculating Lagrangian multipliers. Using the same manipulator properties and taking the inertia tensor information from CAD models, actuator torque values that are required to manipulate upper platform (Eq. (7)) can be drawn in separate graphs (Fig. 13a).

As it can also be seen from Fig. 13b that shows the graphs of differences between the actuator torque values derived by Lagrange method and SimMechanics model, results are consistent with each other. In order to clarify these results, root mean squares (RMS) of actuator torque values derived by Lagrange method and their differences between the SimMechanics model were also compared with each other. At the end promising results of 0.29% and 0.22% deviation was observed (Table IV).

Table V. Dynamixel L54-50-S500-R robot actuator.

Figure 18. Integrated position control diagram taken from the manufacturer documents.

5. Consideration of friction forces

As it can easily be seen from the kinematic structure of the manipulator, upper and lower platforms are sliding through the circular arcs during the operation. This condition will cause friction forces to be generated due to the contacting surfaces. Although these friction forces can be reduced and neglected by utilizing circular bearings or materials with lower friction constants on the real prototype, if necessary they can be implemented to the dynamic analysis of the system. Due to the fact that platform velocity will be low during surgical procedures, inclusion of the friction forces due to the normal forces would be sufficient in the calculations. In light of this, friction forces can be implemented to the Lagrangian formulation via the generalized force component ( $\hat{Q}_{j}$ ) due to an external force. Utilizing rotation matrices of the upper and lower platforms of the manipulator (Eq. (17)), normal components of platform weights can be calculated by projecting them to the platforms unit normal axis as represented by the third columns of rotation matrices.

(20) \begin{equation} \,F_{Nu}=\frac{B_{z}}{r_{u}}m_{u}g,\,\,\,F_{NL}=\frac{B_{z}}{r_{u}}m_{L}g \end{equation}

where in Eq. (20), $F_{Nu}$ represents the normal force component of the upper platforms weight and $F_{Lu}$ represents normal force component of the lower platforms weight. In order to properly find the directions of these friction forces, Eq. (17) should be utilized. It is clear that direction of the friction forces that are generated on individual platforms should be tangent to the related circular rails and opposite to the direction of motion. Unit vectors of these directions are represented by the first column of $\,{}_{u}^{0}{\mathbf{R}}{}$ rotation matrix for the upper platform and by the second column of the $\,{}_{L}^{0}{\mathbf{R}}{}$ rotation matrix for the lower platform. If these are utilized in combination with the signs of platform velocity components in related directions, friction forces due to the normal forces can be expressed as below.

(21) \begin{equation} \,\mathbf{F}_{fu}=\mathrm{sign}\left(-\left[\begin{array}{l@{\quad}l@{\quad}l} \dot{\mathrm{B}}_{x} & \dot{\mathrm{B}}_{y} & \dot{\mathrm{B}}_{z} \end{array}\right]\left[\begin{array}{l} \dfrac{C_{1}\mathrm{B}_{z}-S_{1}\mathrm{B}_{y}}{r_{u}}\\[14pt] \dfrac{S_{1}\mathrm{B}_{x}}{r_{u}}\\[14pt] \dfrac{-C_{1}\mathrm{B}_{x}}{r_{u}} \end{array}\right]\right)\dfrac{\mathrm{B}_{z}}{r_{u}}m_{u}g\mu \left(\dfrac{C_{1}\mathrm{B}_{z}-S_{1}\mathrm{B}_{y}}{r_{u}}\mathbf{i}+\dfrac{S_{1}\mathrm{B}_{x}}{r_{u}}\mathbf{j}+\dfrac{-C_{1}\mathrm{B}_{x}}{r_{u}}\mathbf{k}\right) \end{equation}

(22) \begin{equation} \mathbf{F}_{fL}=\mathrm{sign}\left(-\left[\begin{array}{l@{\quad}l@{\quad}l} \dot{\mathrm{B}}_{x} & \dot{\mathrm{B}}_{y} & \dot{\mathrm{B}}_{z} \end{array}\right]\left[\begin{array}{l} \dfrac{-S_{2}\mathrm{B}_{y}}{r_{u}}\\[14pt] \dfrac{C_{2}\mathrm{B}_{z}+S_{2}\mathrm{B}_{x}}{r_{u}}\\[15pt] \dfrac{-C_{2}\mathrm{B}_{y}}{r_{u}} \end{array}\right]\right)\dfrac{\mathrm{B}_{z}}{r_{u}}m_{L}g\mu \left(\dfrac{-S_{2}\mathrm{B}_{y}}{r_{u}}\mathbf{i}+\dfrac{C_{2}\mathrm{B}_{z}+S_{2}\mathrm{B}_{x}}{r_{u}}\mathbf{j}+\dfrac{-C_{2}\mathrm{B}_{y}}{r_{u}}\mathbf{k}\right) \end{equation}

where $\mu$ represents coefficient of kinetic friction. Substituting Eqs. (21) and (22) into Eq. (7) as generalized force component ( $\hat{Q}_{j}$ ) due to an external force, and selecting proper coefficient of kinetic friction ( $\mu$ ) with respect to the contact surface properties, friction can be implemented to the theoretical dynamic analysis. Friction force values $\mathbf{F}_{fu}$ and $\mathbf{F}_{fl}$ with $\mu =1$ and velocity components on upper and lower rails for the given trajectory can be seen below.

Figure 19. Comparisons of actuator torques and trajectories between virtual and real model (torque 1).

Figure 20. Comparisons of actuator torques and trajectories between virtual and real model (torque 2).

Figure 21. Dual motion capture cameras and manipulator verification setup.

Figure 22. Attached reference markers on the manipulator with landmark positions.

Figure 23. Stylus positioning on landmark positions during measurement procedure.

Figure 24. IR reflector attached on the platform and motion execution.

It should be noted that in order to acquire a clear representation, coefficient of kinetic friction ( $\mu$ ) was taken as 1 in Fig. 14. This value will be 0.002 to 0.003 in case of utilization of circular bearings, 0.04 in case of Teflon on Teflon contact and 0.24 in case of Teflon on aluminum contact.

Figure 25. Comparisons between actual and desired trajectories.

6. Preliminary hardware verification

In order to verify reliability of the proposed virtual model and dynamic analysis, similar but faster trajectory (Fig. 15) of the case study was tracked by the manufactured version of the updated spherical manipulator (Fig. 16).

As shown in Fig. 16, main parts of the manipulator prototype (actuator shafts, circular links, base) were manufactured by using aluminum and housings were manufactured via rapid prototyping by using ABS filament in order to reduce overall weight of the system without losing structural integrity. As proposed systems platforms slide on circular links, in order to reduce overall friction as much as possible during the operation, Teflon guides were manufactured and assembled (Fig. 17).

After all of the necessary parameters were updated in virtual environment considering manufactured prototype of the manipulator, new joint trajectories (Fig. 15) were tracked by the virtual manipulator model and real prototype that was actuated by two Dynamixel L54-50-S500-R robot actuators simultaneously in order to acquire torque values. Physical properties and dimensions of the actuators can be seen in Table V.

Figure 26. Centroids of the two data sets.

Friction force between Teflon guides and aluminum rails was simulated in the virtual model by increasing joint damping constants in SimMechanics. In order to control robot manipulator throughout the study, only the integrated position controller of the actuators provided by the manufacturer was utilized. Desired joint trajectories (Fig. 15) that were calculated for the hardware verification trials were fed to the actuator position controller as goal positions. Integrated position control diagram taken from the manufacturer (Robotis) documents can be seen in Fig. 18.

After the trials, results that were extracted from the virtual model and L54-50-S500-R actuators were compared for verification purposes. Due to the fact that actuators used in the system have no actual torque sensor, torque readings were calculated by acquiring current values during the operation and mapping them to the current-torque actuator performance graph given by the manufacturer.

As seen in Figs. 19 and 20, although acquired current data received from the actuators are noisy, results that were extracted from the virtual model and L54-50-S500-R actuators are consistent with each other. In light of this, it can be concluded that constructed virtual model in MATLAB SimMechanics for the actual hardware gives promising results in terms of actuator torque values.

In order to strengthen the case and proposed application feasibility, positioning precision of the manipulator was also verified by using motion capture cameras in order to be sure that desired trajectories were actually tracked by the manipulator platform with minimal error. For this task, dual OptiTrack V100R2 motion capture cameras were utilized (Fig. 21).

Prior to the experiments, in order to find relationship between motion capture setup measurement reference (C) frame with manipulator reference frame (M) by calculating transformation matrix ${}_{C}^{M}{\mathbf{T}}{}$ [Reference Uslu, Guzin and Gezgin22,Reference Uslu, Gezgin, Özbek, Güzin, Can and Çetin23], small reference landmark markers (Fig. 22) were designed, manufactured, and attached to the manipulator with precisely known landmark positions ${}^{M}{\!\mathbf{\rho }}{_{i}^{}}(i=1,\ldots,8)$ .

Following the camera calibration and the registration procedure given in detail in ref. [Reference Uslu, Guzin and Gezgin22–Reference Uslu, Gezgin, Özbek, Güzin, Can and Çetin23], measurements from landmark positions were captured by using calibrated stylus with attached IR reflectors (Fig. 23).

Using measured position data ${}^{C}{\!\mathbf{\rho }}{_{i}^{}}(i=1,\ldots,8)$ and actual landmark positions ${}_{}^{M}{\mathbf{\rho }}{_{i}^{}}(i=1,\ldots,8)$ transformation matrix ${}_{C}^{M}{\mathbf{T}}{}$ was calculated with an RMS error of 0.587 mm [Reference Uslu, Guzin and Gezgin22–Reference Uslu, Gezgin, Özbek, Güzin, Can and Çetin23]. In order to compare desired end effector trajectory (Fig. 15) with actual trajectory, an IR reflector was attached to the center of systems platform (Fig. 24), and its path was captured by the calibrated motion capture setup during the motion.

Once extracted data were transformed to the manipulator reference via the transformation matrix ${}_{C}^{M}{\mathbf{T}}{}$ calculated before, desired and actual trajectories of the manipulator were compared with each other (Fig. 25).

As seen in Fig. 25, trajectory tracking performance and positioning precision of the manufactured prototype of the manipulator give promising results. Between the centroids of the tracked data and desired trajectory data, 0.731 mm deviation was calculated (Fig. 26).

Considering the fact that actual RMS error of the registration was 0.587 mm calculated deviation was meaningful.