Development of a spherical parallel manipulator for brain surgery applications: preliminary study on the dynamic analysis and verification

Abstract

This study focuses on development, task planning, and dynamic analysis of a previously proposed spherical parallel robot manipulator that is conceptually enhanced to adapt various brain surgery scenarios. Conceptual design of the proposed manipulator was briefly introduced and explained. In order to simulate one of the possible surgery scenarios, a case study of craniotomy was designed along with its trajectory planning. Dynamic analysis of proposed manipulator was performed by Lagrange method, and required actuator torque values were calculated for the desired trajectory. At the end of the study, hardware verification was carried out on the manufactured prototype of the system by comparing both calculated/acquired torque values and desired/actual trajectories. Promising verification results in terms of system dynamics and trajectory execution were introduced.

1. Introduction

Advanced robot technologies have been developed rapidly throughout the world in conjunction with the utilization of robot assistance in medical field. Due to its potentials to improve precision, enhance dexterity, eliminate infection risks, and reduce complication rates, robot- assisted surgeries have continually received wide attention in both academic literature and medical industry. Over the past decade, various applications of medical robotics were proposed and implemented to the field of medicine in order to improve traditional operations by the assistance of robotic systems. In order to achieve required effectiveness to fulfill delicate surgical tasks, these robotic applications mainly depend on efficient operation planning and verification via different methodologies. Prior to the initiation of surgical operation trials, verifications can be carried out by utilizing various modeling techniques and simulation tools. During these processes, end effector trajectories and system behaviors can be observed in a virtual environment and at the same time they can be improved with respect to the acquired information. This fact not only helps further development of manipulators by providing knowledge of system characteristics but also helps surgeon by means of enhanced control over the actual system via advanced visualization, targeting, and task execution.

In order to acquire system behavior throughout operation planning or task formation, end effector trajectories are generated by utilizing consecutive via points that are located inside the workspace of the manipulator and tracked by actual system or its virtual model. Thus, derivation of kinematic and dynamic modeling along with simulation visualization plays an important role to acquire desired information. Gosselin et al. [1] proposed a dynamic trajectory planning method for three degrees of freedom cable suspended parallel robot manipulator. Kinematic and dynamic modeling of three cable spatial manipulator that has a point mass end effector was carried out. The authors introduced a technique for the generation of point to point trajectories, where consecutive points that form trajectory segments defined by trigonometric functions have zero instantaneous velocity and ensure acceleration continuity. Proposed method also includes application of natural frequencies and was verified on a manufactured manipulator with an example trajectory. Zha [2] introduced an approach for manipulator pose trajectory planning by utilizing genetic algorithms and compared proposed approach with several classical optimization algorithms in a MATLAB environment. Case studies were carried out and evaluated on a 3R planar and PUMA 560 serial robot manipulators by considering their kinematic and dynamic constraints. Acquisition of system characteristics or behaviors during end effector trajectory tracking becomes superior by the derived model of manipulator kinematics and dynamics [3]. Several approaches [4] are available to characterize dynamics of robot manipulators as virtual work principle [5], Newton–Euler and Lagrange’s formulations. Due to the fact that Lagrange formulation eliminates the consideration of reaction forces and moments, it reduces procedural complexity to derive manipulators equation of motion and thus preferred by many studies. Elkady et al. [6] and Yen et al. [7] studied on the dynamic modeling of three degrees of freedom Cartesian manipulator. Both of the groups derived dynamics of the manipulator by utilizing Lagrangian methods. Various types of motion controllers were tried on a MATLAB simulation environment by utilizing a given end effector trajectory and acquired performances were compared throughout the studies. Guo et al. [8] presented the dynamic analysis of the Gough-Stewart manipulator by using a methodology that includes utilization of both Newton–Euler and Lagrange formulation along with manipulator kinematics. The authors verified the effectiveness of the proposed method via simulations in MATLAB software. Cheng et al. [9] proposed a methodology for the derivation of dynamics for both normal and redundant parallel manipulators via Lagrange-D’Alembert formulation. Experimental simulation results were given by the authors for various control methods throughout their study. Considering studies of both trajectory tracking and derivation of dynamic characteristics of any manipulator, it can easily be seen that there exist a necessity of a simulation medium for methodological verification and performance evaluation. Throughout the related literature [10–13], various types of mediums have been utilized for these purposes. Adam et al. [14] carried out kinematic and dynamic modeling of a 5R robot manipulator. The authors verified their results in a virtual CATIA and MATLAB simulation environment by using a formed trajectory. Joint torque values were compared for both elbow up and elbow down postures of the manipulator. Similarly, Gouasmi et al. [15] compared dynamic behavior of a 2R serial manipulator under a given trajectory for both manipulator postures on MATLAB simulink and solidworks simulations to verify the reliability of proposed models.

In light of related literature and by also considering the studies of robot assistance in brain surgery operations [16–18], this paper first presents further development of two degrees of freedom spherical parallel manipulator that was first introduced by the study of same authors [19]. Throughout the study, enhanced modular conceptual design of a manipulator to adapt various brain surgery operations was briefly introduced. In order to simulate manipulator dynamics under given task, virtual model of the system was constructed on MATLAB SimMechanics environment. A case study that involves surgical application of opening a bone flap on human skull was selected. Trajectory planning of the selected case was performed and necessary actuator torque values were extracted from the virtual model of the manipulator. In order to verify reliability of the model, dynamic analysis of the manipulator was carried out by using Lagrange method and detailed information about the procedure was introduced. Results of the virtual simulation and theoretical approach were compared for verification purposes. Last part of the study was dedicated to preliminary hardware verification. Desired trajectory of the case study was tracked by the manufactured prototype of the proposed system to extract actuator torque requirements and precision performance of the system. Acquired results were compared and discussed at the end of the paper.

2. Structural design

As mentioned in previous study of the authors [19], 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 [19], 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 [20] 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 [21] 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) [21].

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. [19] (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}

3.2. SimMechanics simulation

After calculating actuator joint angle functions, they were applied to the imported manipulator model inside SimMechanics (Fig. 9a), and torque values of each actuators were extracted from the environment during simulation run (Fig. 9b). Each torque graph depicts necessary actuator torque that is required to manipulate upper platform of the system with respect to the given trajectory.

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 [4], 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}}{}$ [22,23], 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. [22–23], 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 [22–23]. 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.

7. Conclusion

This study introduced a conceptual design of brain surgery manipulator that has been enhanced from two degrees of freedom spherical parallel manipulator presented in the earlier study of same authors [19]. Addition of different modules to the platform of parallel manipulator has allowed possible adaptation of the proposed system to various brain surgery operations such as craniotomy, neuroendoscopy, and deep brain stimulation. In order to demonstrate one of the mentioned possible tasks, a case study of opening bone flap on human skull was selected and trajectory planning of the surgical operation was introduced. Virtual model of the manipulator was constructed in MATLAB SimMechanics environment to reveal dynamic behavior of the system under given trajectory by means of actuator torque requirements. Reliability of the virtual model was verified by carrying out manipulators analytical dynamic analysis by using Lagrange method. Procedural approach of the dynamic analysis was also introduced in detail. Comparison between virtual and analytical model revealed a promising 0.29% and 0.22% deviation for the required torque values of actuators, respectively. At the end of the study, constructed virtual model and proposed dynamic analysis procedure were also confirmed via utilizing hardware verification procedure on a manufactured prototype of the spherical manipulator along with its trajectory tracking performance.

Considering selected case study and comparisons between virtual and Lagrange model, promising results were acquired throughout the study that validates the applicability of both methodologies in a dedicated brain surgery application. In light of this, current paper contributes related literature not only by proposing a modified conceptual design of a brain surgery manipulator but also proposing a methodological approach for task planning along with the validation of manipulator dynamics and precision.