**New Frontiers in Parallel Robots**

Editors

**Zhufeng Shao Dan Zhang St ´ephane Caro**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Zhufeng Shao Tsinghua University China

Dan Zhang York University Canada

Stephane Caro ´ National Centre for Scientific Research (CNRS) and the Laboratory of Digital Sciences of Nantes (LS2N) France

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Machines* (ISSN 2075-1702) (available at: https://www.mdpi.com/journal/machines/special issues/Parallel Robots).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-7252-9 (Hbk) ISBN 978-3-0365-7253-6 (PDF)**

Cover image courtesy of Dan Zhang and Stephane Caro ´

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Editors**

#### **Zhufeng Shao**

Zhufeng Shao, Ph.D., is an associate professor in the Department of Mechanical Engineering at Tsinghua University, a member of the National Industrial Foundation Expert Committee, a special expert of the Chinese Academy of Engineering for foundation strategy research, a senior member of the Chinese Mechanical Engineering Society, an Editorial Board member of the journal *Defence Technology* and a member of the National Professional Standards Technical Committee. He has won the Silver Award of China Patent, the First Prize of Technical Invention of the Science and Technology Award of the China Instrument Society, the Special Prize and Second Prize of Scientific and Technological Progress of the China Machinery Industry Science and Technology Award and the Second Prize of Natural Science of the Ministry of Education, as well as others. Dr. Shao received his Ph.D. in Mechanical Engineering from Tsinghua University, China, in June 2011. He is mainly engaged in the research of high-performance robotic equipment and intelligent manufacturing. Dr. Shao's research interests include parallel mechanisms, cable-driven parallel robots, soft robots, the innovative design of robots, macro and micro robot systems, humanoid and rehabilitation robots, autonomous systems and intelligent manufacturing. He has published 38 SCI-indexed papers, 32 EI-indexed papers and two books, and has been granted 45 Chinese invention patents and 12 software copyrights. He has developed an intelligent painting system for large civil aircraft and an automatic painting system for the non-structural surface of large ship segments, and has invented a new type of cable-driven high-speed parallel robot. Dr. Shao led the construction of a remote operation and maintenance system for CNC metal-cutting machine tools and has released six national standards related to intelligent manufacturing.

#### **Dan Zhang**

Dan Zhang, Ph.D., is a Kaneff Professor and Tier 1 York Research Chair in Advanced Robotics and Mechatronics in the Department of Mechanical Engineering at York University. Dr. Zhang was a Canada Research Chair in Advanced Robotics and Automation and was a founding Chair of the Department of Automotive, Mechanical and Manufacturing Engineering with the Faculty of Engineering and Applied Science at Ontario Tech University. He received his Ph.D. in Mechanical Engineering from Laval University, Canada, in June 2000. Dr. Zhang's research interests include the synthesis and optimization of parallel and hybrid mechanisms, generalized parallel mechanisms research, reconfigurable robots, the innovative design of parallel robots, the parallelization of serial robots, micro/nano manipulation and mems devices (e.g., sensors), rescue robots, smart biomedical instruments (e.g., exoskeleton robots and rehabilitation robotics), AI/robotics/autonomous systems, aerial and underwater robotics, artificial intelligence for robotics and intelligent reconfigurable adaptive landing gear and manipulators (manipulanders). Dr. Zhang has published 241 journal papers and 187 conference papers, 12 books, nine book chapters and numerous other technical publications. Dr. Zhang has served as a General Chair for 67 International Conferences and delivered 117 keynote speeches. Dr. Zhang was listed in the World's Top Two Percent Researchers by Stanford's Standardized Citation Indicators in 2020 and 2021. Dr. Zhang is a Fellow of the Canadian Academy of Engineering (CAE), a Fellow of the Engineering Institute of Canada (EIC), a Fellow of the American Society of Mechanical Engineers (ASME), a Fellow of the Canadian Society for Mechanical Engineering (CSME), a Senior Member of the Institute of Electrical and Electronics Engineers (IEEEs) and a Senior Member of SME.

#### **St ´ephane Caro**

Stephane Caro, Ph.D., is the Director of Research at the National Centre for Scientific Research ´ (CNRS) and at the Laboratory of Digital Sciences of Nantes (LS2N), UMR CNRS 6004, France. Dr. Caro received his Engineering and M.Sc. degrees in mechanical engineering from Ecole Centrale Nantes, Nantes, France, in 2001, and his Doctorate degree in mechanical engineering from the University of Nantes in 2004. He was a Post-Doctoral Fellow in the Centre for Intelligent Machines, McGill University, Montreal, QC, Canada, from 2005 to 2006. He was awarded the accreditation to supervise research (HDR) in 2014. Dr. Caro's research focuses on the design, modeling and control of cable-driven parallel robots and reconfigurable parallel robots. Dr. Caro is the head of the "Robots and Machines for Manufacturing Society and Services" (RoMaS) Team at LS2N. He is also a part-time researcher at IRT Jules Verne, a mutualized industrial research institute, and a lecturer at Ecole Centrale de Nantes. Dr. Caro is a member of ASME and IEEE and was one of the two recipients of the 2018 Reviewers of the Year for the ASME *Journal of Mechanisms and Robotics*. He was also one of two recipients of the 2019 Crossley Award for Mechanism and Machine Theory. Dr. Caro is Associate Editor for the ASME *Journal of Mechanisms and Robotics* and the *IEEE Robotics and Automation Letters*. He is a member of the IFToMM Computational Kinematics and ASME Mechanisms and Robotics Technical Committees.

### *Editorial* **New Frontiers in Parallel Robots**

**Zhufeng Shao 1,\*, Dan Zhang 2 and Stéphane Caro 3**

In the field of parallel robots, marked by the birth and application of the Gough– Stewart parallel mechanism [1] in the 1960s, great progress has been made in the past 60 years. The most notable feature of a parallel robot is that there are multiple closed-loop branch chains jointly connecting and driving the moving platform [2], which gives great flexibility in its configurations, creating a new way to change performance through robot configuration. Parallel robots usually have outstanding advantages of high stiffness, high precision, and high speed [3], which make up for the performance shortcomings of serial robots. The abundant configurations and complex mechanisms of parallel robots also present challenges in terms of configuration synthesis, performance evaluation, modeling, calibration, control, etc. Opportunities coexist with challenges, and parallel robots attract attention from both academia and industry. Today, parallel robots are constantly enriched, and new types of parallel robots, such as cable-driven parallel robots (CDPRs) [4], soft parallel robots [5], and hybrid robots [6], are constantly emerging. In particular, while inheriting the abovementioned advantages, a CDPR has the advantages of low cost, high energy efficiency, easy reconfiguration, and light weight, showing great application potential in scenarios such as large working spaces, heavy loads, high speeds, and bionics [7,8]. Parallel robotics research and applications show continued vitality and are expected to transform the industry in the future.

This Special Issue provides an international forum for professionals, academics, and researchers to present the latest developments from theoretical studies and applications of parallel robots. It includes 10 selected papers, covering important aspects of parallel robots such as modeling and control, error analysis and calibration, singularity analysis, and trajectory planning. The contents of these studies are briefly described here.

In [9], an evaluation model is established to analyze the influence of geometric errors on limbs' comprehensive deformations for an over-constrained parallel manipulator. The evaluation model is established based on kinematics and verified through simulations. Two global sensitivity indices are proposed, and a sensitivity analysis is conducted using the Monte Carlo method throughout the reachable workspace. The geometric errors that have greater effects on the average angular comprehensive deformation are identified.

In [10], a consistent solution strategy for static equilibrium workspaces of different types of under-constrained cab-driven robots is presented. The dynamic models and parameters that are applied to make the system stable for point-to-point movements are introduced. The constraints of the dynamics model are incorporated into the trajectory planning process to achieve point-to-point trajectory planning for the under-constrained cable-driven robots.

In [11], the authors present the singularity analysis and the geometric optimization of a 6-DOF (Degrees of Freedom) parallel robot for SILS (Single-Incision Laparoscopic Surgery). Based on a defined set of input/output constraint equations, the singularities of the parallel robotic system are determined and geometrically interpreted. Then, the

**Citation:** Shao, Z.; Zhang, D.; Caro, S. New Frontiers in Parallel Robots. *Machines* **2023**, *11*, 386. https:// doi.org/10.3390/machines11030386

Received: 13 March 2023 Accepted: 13 March 2023 Published: 15 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

geometric parameters for the 6-DOF parallel robot are optimized to make the operational workspace singularity-free.

Paper [12] focuses on pick-and-place trajectory planning and tracking control of a cable-based gangue-sorting robot in the operation space. A four-phase pick-and-place trajectory planning scheme based on an S-shaped acceleration/deceleration algorithm and the quantic polynomial trajectory planning method is proposed. A robust adaptive fuzzy tracking control strategy is presented against inevitable uncertainties and unknown external disturbances. The proposed method guarantees a stable and accurate pick-andplace trajectory tracking process.

Paper [13] proposes a fractional-order impedance control scheme, named KDHD, in which additional damping is added, proportional to the half-order derivatives of the endeffector position errors according to the half-derivative damping matrix, HD. The proposed impedance controller represents an extension to multi-input multi-output robotic systems of the PDD1/2 controller for single-input single-output systems, which over performs the PD scheme in the transient behavior.

Paper [14] focuses on the dynamic modeling, workspace analysis, and multi-objective structural optimization of a large-span, high-speed, cable-driven parallel camera robot. The curved cable, due to the self-weight, is modeled as a catenary, and the dynamic model is derived by decomposing the motion of the cable into an in-plane motion and an out-plane motion. An optimization model is presented to simultaneously improve the workspace volume, anti-wind disturbance ability, and impulse of tensions on the camera and pan-tilt device system (CPTDS).

In [15], the authors present kinematic and dynamic modeling and workspace analysis for a novel suspended CDPR which generates Schönflies motions. The kinematics of the CDPR are solved through a geometrical approach. The dynamic feasible workspace of the robot is determined. Experiments are performed on a prototype of the robot to demonstrate the correctness of the derived models and workspace.

Paper [16] proposes a new method for the kinematic calibration of parallel robots to strict pose error bounds. The new method includes a new pose error model with 60 error parameters and a different kinematic parameter error identification algorithm based on Linfinity parameter estimation. Parameter errors are identified by using linear programming. The feasibility and validity of the proposed kinematic calibration are verified through both simulations and experiments.

Paper [17] studies the 3-DOF cutting stability and surface quality optimization of a parallel kinematic manipulator (PKM). A prediction model for the 3-DOF stability of helical milling based on the PKM is established through a semi-discrete method based on the natural frequency analysis of the PKM and a cutting force model of titanium alloy helical milling. A step-cutter is used to improve the machining process by enhancing the stability domain. The proposed method can provide a reference for further optimization of the prediction and optimization of the milling of difficult-to-process materials based on a PKM in the future.

In [18], the authors develop a simple model to evaluate the first natural frequencies of over-constrained PKMs. The PKM legs are modeled by beams, and constraint equations between the parameters are determined according to screw theory. The focus of this paper is to determine the global mass and stiffness matrices of the PKM in stationary configurations without the use of Jacobian matrices. The proposed method can be easily used at the conceptual design stage of PKMs.

The Guest Editors thank all of our colleagues who have taken interest in this Special Issue, especially the authors of the papers published in this Special Issue. All the of papers underwent a rigorous review process to ensure the high quality of the publications. We are grateful to the reviewers who evaluated these papers and provided valuable comments based on their professional perspectives. We also would like to thank the editors from MDPI for their support and effort in the organization and publication of this Special Issue.

It is hoped that the papers published in this Special Issue can be used as vehicles to promote knowledge sharing in the field of parallel robots. More importantly, we hope more people will be informed about and understand parallel robots and their latest technologies and actively participate in the innovation, research, development, and application promotion of parallel robots.

**Author Contributions:** Conceptualization, Z.S., D.Z. and S.C.; writing—original draft preparation, Z.S.; writing—review and editing, D.Z. and S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. U19A20101) and the National High-tech Ship Research Project of China (No. MC-202003-Z01).

**Acknowledgments:** We express great thanks to the editors of MDPI for their excellent support for this Special Issue, and it would have been impossible without their persistence and help. Many thanks to all our colleagues who are interested in these research topics and submitted their research for our Special Issue. A special thanks to the reviewers for their efforts and time spent in order to maintain the high quality of all contributions.

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

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Geometric Error Analysis of a 2UPR-RPU Over-Constrained Parallel Manipulator**

**Xu Du 1, Bin Wang 1 and Junqiang Zheng 2,\***

**Abstract:** For a 2UPR-RPU over-constrained parallel manipulator, some geometric errors result in internal forces and deformations, which limit the improvement of the pose accuracy of the moving platform and shorten the service life of the manipulator. Analysis of these geometric errors is important for restricting them. In this study, an evaluation model is established to analyse the influence of geometric errors on the limbs' comprehensive deformations for this manipulator. Firstly, the nominal inverse and actual forward kinematics are analysed according to the vector theory and the local product of the exponential formula. Secondly, the evaluation model of the limbs' comprehensive deformations is established based on kinematics. Thirdly, 41 geometric errors causing internal forces and deformations are identified and the results are verified through simulations based on the evaluation model. Next, two global sensitivity indices are proposed and a sensitivity analysis is conducted using the Monte Carlo method throughout the reachable workspace of the manipulator. The results of the sensitivity analysis indicate that 10 geometric errors have no effects on the average angular comprehensive deformation and that the identified geometric errors have greater effects on the average linear comprehensive deformation. Therefore, the distribution of the global sensitivity index of the average linear comprehensive deformation is more meaningful for accuracy synthesis. Finally, simulations are performed to verify the results of sensitivity analysis.

**Keywords:** 2UPR-RPU parallel manipulator; over-constrained parallel manipulator; geometric error; deformation; sensitivity analysis

#### **1. Introduction**

Parallel mechanisms with three DOFs have been successfully applied to hybrid serial– parallel machine tools, such as the well-known Eco-speed series, Tricept, and Exechon [1–6], owing to their high stiffness, large payload, and good dynamics. To achieve a simpler structure, Li et al. [3] designed a 2R1T (R denotes a rotational DOF, and T denotes a translational DOF) parallel mechanism named 2UPR-RPU. This mechanism is not only easier to control but also suitable for many operations along the surfaces. However, it is an over-constrained parallel mechanism with common constraints and over-constraints [7,8]. Some geometric errors in a manipulator based on this mechanism break the common constraints and over-constraints, resulting in internal forces and deformations. The internal forces and deformations not only limit the further improvement of the pose accuracy of the moving platform but also shorten the service life of the manipulator [9,10]. Therefore, it is necessary to restrict the internal-force-and-deformation-related geometric errors in the 2UPR-RPU parallel manipulator.

The accuracy design [11–13] can be applied to restrict geometric errors by determining the tolerances of the fabrication and assembly of machines. It consists of three components: error modelling [14–16], sensitivity analysis [17–19], and accuracy synthesis [20–22], where error modelling is the basis of sensitivity analysis and accuracy synthesis. Zhang et al. [13] applied the closed-loop vector and first-order perturbation methods to establish a geometric

**Citation:** Du, X.; Wang, B.; Zheng, J. Geometric Error Analysis of a 2UPR-RPU Over-Constrained Parallel Manipulator. *Machines* **2022**, *10*, 990. https://doi.org/10.3390/ machines10110990

Academic Editors: Zhufeng Shao, Dan Zhang and Stéphane Caro

Received: 18 September 2022 Accepted: 27 October 2022 Published: 29 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

error model for a 2UPR-RPS over-constrained manipulator, and they identified the geometric errors that affected the pose errors of the moving platform. Zhang et al. [15] utilised the screw theory to establish a geometric error model for a 4RSR-SS over-constrained parallel tracking machine. With the use of the geometric error model, 53 geometric errors that had a significant influence on the pose errors of the moving platform were identified after sensitivity analysis. However, neither of the above two methods considers the deformations caused by internal forces in over-constrained parallel manipulators. Taking parameter uncertainties into account, Tang et al. [23] built a general interval kinetostatic model for a 2UPR-SPR over-constrained parallel machine to perform sensitivity analysis and tolerance allocation. To predict the pose errors of an over-constrained extendible support structure, Yu et al. [24] proposed a comprehensive model that simultaneously considered geometric errors, joint gaps, and link flexibility. In spite of good accuracy, these two models are complicated for the stiffness matrix needs to be derived and the stiffness coefficients of parts need to be obtained via finite element software.

Affected by geometric errors, the end poses of different limbs of a parallel manipulator should be theoretically inconsistent. However, they can be consistent in nonoverconstrained parallel manipulators due to the existence of the moving platform and the motion deviations of passive joints. On this basis, a numerical iterative algorithm [25,26] was proposed to analyse the kinematics of non-overconstrained parallel manipulators with kinematic errors. Inspired by this algorithm, this study aims to establish an evaluation model based on kinematics to analyse the influence of geometric errors on the limbs' comprehensive deformations for the 2UPR-RPU over-constrained parallel manipulator.

Based on the established evaluation model, sensitivity analysis can help reveal the influence of different internal-force-and-deformation-related geometric errors on the limbs' comprehensive deformations. The interval analysis method and probabilistic method have been commonly used for sensitivity analysis of the moving platform's pose error in literature. The interval analysis method treats geometric errors as interval variables and can get a balance between calculation speed and accuracy [11,18]. Treating geometric errors as random variables with a normal distribution, the probabilistic method can be divided into the Monte Carlo method and the probability modelling method. The Monte Carlo method calculates the moving platform's pose errors according to the geometric error model and lots of random values of a geometric error [22,27]. It has good accuracy and low computational efficiency. The probability modelling method establishes an analytical model between the standard deviation of each geometric error and that of the moving platform's pose error based on the geometric error model [28]. In spite of high computational efficiency, this method needs prior knowledge about probability distributions. Considering that the interval analysis method and probability modelling method are not suitable when the geometric error model is iterative, the Monte Carlo method is utilised to analyse the influence of geometric errors on the limbs' comprehensive deformations in this paper.

The remainder of this paper is organised as follows. In Section 2, the 2UPR-RPU parallel mechanism is briefly introduced. Section 3 presents an analysis of the nominal inverse kinematics and actual forward kinematics. Section 4 establishes an evaluation model of the limbs' comprehensive deformations caused by geometric errors. Based on the evaluation model, the internal-force-and-deformation-related geometric errors are identified and the results are verified through simulations in Section 5. In Section 6, two global sensitivity indices are proposed and sensitivity analysis is conducted. Simulations are also performed to verify the results of sensitivity analysis. Finally, the conclusions are drawn in Section 7.

#### **2. 2UPR-RPU Parallel Mechanism**

As shown in Figure 1, the 2UPR-RPU parallel mechanism mainly consists of a moving platform, two UPR limbs, one RPU limb, and one fixed base, where the moving platform and fixed base are represented by the isosceles right triangles ΔA1A2A3 and ΔB1B2B3. U, P, and R denote universal, prismatic, and revolute joints, respectively. **B**1, **B**2 and **A**3 are the centres of U, and **A**1, **A**2 and **B**3 are the centres of R. Because each universal joint is

equivalent to two mutually perpendicular revolute joints, the UPR limb is equivalent to the RRPR limb, and the RPU limb is equivalent to the RPRR limb. The axis of the *j*th joint of the *i*th limb is denoted by **s***i*,*j*. A fixed coordinate system {**o**B; **x**, **y**, **z**} is established at the midpoint between **B**1 and **B**2, where **x** points from **B**2 to **B**1 and **y** points from **o**B to **B**3. Similarly, a moving coordinate system {**o**A; **u**, **v**, **w**} is also established, where **u** points from **A**2 to **A**1 and **v** points from **o**A to **A**3. The coordinate axes **z** and **w** are determined using the right-hand rule. For the 2UPR-RPU parallel mechanism, each limb exerts a force and a couple on the moving platform [8], where the two forces from the UPR limbs are parallel to **v**, and the three couples from the UPR and RPU limbs rotate around **w**. It is worth mentioning that the two forces parallel to **v** will lead to over-constraint, and the three couples rotating around **w** will lead to common constraints. Thus, the 2UPR-RPU parallel mechanism is an over-constrained parallel mechanism.

**Figure 1.** Schematic diagram of the 2UPR-RPU parallel mechanism.

#### **3. Kinematics**

Inverse kinematics aims to calculate the displacements of all joints relative to their initial positions or angles according to a given target pose of the moving platform. Forward kinematics is the reverse operation of inverse kinematics. Inverse kinematics without considering geometric errors is called nominal inverse kinematics. In this section, the nominal inverse kinematics of actuated joints and passive joints is first introduced. Then, the actual forward kinematics of the limbs is derived.

#### *3.1. Nominal Inverse Kinematics*

The position and orientation of the moving platform shown in Figure 1 can be described by - *xyz*T and - *αβγ*T , respectively, where - *xyz*T denotes the position coordinates of **o**A with respect to {**o**B; **x**, **y**, **z**} and - *αβγ*T denotes the Euler angle with respect to **z-x**-**v**. Because only the translation motion along **o**B**o**A and the rotations around **x** and **v** can be achieved by the moving platform [8], - *z β γ*T is sufficient to represent

the poses. For a given target pose of the moving platform, the nominal displacements of actuated P-joints can be derived using the closed-loop vector method [10] as follows:

$$\begin{cases} \begin{aligned} q\_{1,3} &= ||\mathbf{B}\_1 \mathbf{A}\_1|| - ||\mathbf{B}\_1 \tilde{\mathbf{A}}\_1|| \\ q\_{2,3} &= ||\mathbf{B}\_2 \mathbf{A}\_2|| - ||\mathbf{B}\_2 \tilde{\mathbf{A}}\_2|| \\ q\_{3,2} &= ||\mathbf{B}\_3 \mathbf{A}\_3|| - ||\mathbf{B}\_3 \tilde{\mathbf{A}}\_3|| \end{aligned} \end{cases} \tag{1}$$

where -· represents the Euclidean norm. **~ A***i* denotes the initial position of **A***i*, which is determined by

$$\begin{cases} \begin{array}{ll} \mathbf{B}\_{1}\mathbf{A}\_{1} = \begin{bmatrix} \ l\_{\text{A}}\cos\gamma - l\_{\text{B}} & l\_{\text{A}}\sin\beta\sin\gamma - z\tan\beta & -l\_{\text{A}}\cos\beta\sin\gamma + z \end{bmatrix}^{\text{T}}\\ \mathbf{B}\_{2}\mathbf{A}\_{2} = \begin{bmatrix} - & l\_{\text{A}}\cos\gamma + l\_{\text{B}} & -l\_{\text{A}}\sin\beta\sin\gamma - z\tan\beta & l\_{\text{A}}\cos\beta\sin\gamma + z \end{bmatrix}^{\text{T}}\\ \mathbf{B}\_{3}\mathbf{A}\_{3} = \begin{bmatrix} 0 & l\_{\text{A}}\cos\beta - z\tan\beta - l\_{\text{B}} & l\_{\text{A}}\sin\beta + z \end{bmatrix}^{\text{T}} \end{array} \tag{2}$$

where *l*A = -**A**1**A**2-/2 and *l*B = -**B**1**B**2-/2.

For the first UPR limb, the first, second, and fourth joints are passive. The nominal displacement of the first joint can be expressed as

$$q\_{1,1} = \beta \tag{3}$$

The nominal displacement of the second joint can be expressed as

$$q\_{1,2} = \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_1 \mathbf{A}\_1}{||\mathbf{B}\_1 \mathbf{A}\_1||}\right) - \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_1 \tilde{\mathbf{A}}\_1}{||\mathbf{B}\_1 \tilde{\mathbf{A}}\_1||}\right) \tag{4}$$

where **e**1 is the unit vector along **x**.

The nominal displacement of the fourth joint can be expressed as

$$q\_{1,4} = \arccos\left(\frac{-\left(\mathbf{A}\_1\tilde{\mathbf{A}}\_2\right)^T \left(\mathbf{B}\_1\tilde{\mathbf{A}}\_1\right)}{||\mathbf{A}\_1\tilde{\mathbf{A}}\_2|| ||\mathbf{B}\_1\tilde{\mathbf{A}}\_1||}\right) - \arccos\left(\frac{-\left(\mathbf{A}\_1\mathbf{A}\_2\right)^T \left(\mathbf{B}\_1\mathbf{A}\_1\right)}{||\mathbf{A}\_1\mathbf{A}\_2|| ||\mathbf{B}\_1\mathbf{A}\_1||}\right) \tag{5}$$

Because the two UPR limbs are symmetrically distributed with respect to **o**A**o**B, we have

$$q\_{2,1} = \beta \tag{6}$$

$$q\_{2,2} = \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_2 \mathbf{A}\_2}{||\mathbf{B}\_2 \mathbf{A}\_2||}\right) - \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_2 \tilde{\mathbf{A}}\_2}{||\mathbf{B}\_2 \tilde{\mathbf{A}}\_2||}\right) \tag{7}$$

$$q\_{2,4} = \arccos\left(\frac{(\mathbf{A}\_1\mathbf{A}\_2)^T(\mathbf{B}\_2\mathbf{A}\_2)}{||\mathbf{A}\_1\mathbf{A}\_2|| \, ||\mathbf{B}\_2\mathbf{A}\_2||}\right) - \arccos\left(\frac{\left(\mathbf{A}\_1\tilde{\mathbf{A}}\_2\right)^T\left(\mathbf{B}\_2\tilde{\mathbf{A}}\_2\right)}{||\mathbf{A}\_1\tilde{\mathbf{A}}\_2|| \, ||\mathbf{B}\_2\tilde{\mathbf{A}}\_2||}\right) \tag{8}$$

Similarly, the nominal displacements of the first, third, and fourth joints of the RPU limb can be expressed as

$$q\_{3,1} = \arccos\left(\frac{-\mathbf{e}\_2^T \mathbf{B}\_3 \mathbf{A}\_3}{||\mathbf{B}\_3 \mathbf{A}\_3||}\right) - \arccos\left(\frac{-\mathbf{e}\_2^T \mathbf{B}\_3 \tilde{\mathbf{A}}\_3}{||\mathbf{B}\_3 \tilde{\mathbf{A}}\_3||}\right) \tag{9}$$

$$q\_{3,3} = \arccos\left(\frac{(\mathbf{Re}\_2)^T(\mathbf{B}\_3\mathbf{A}\_3)}{||\mathbf{B}\_3\mathbf{A}\_3||}\right) - \arccos\left(\frac{\left(\tilde{\mathbf{Re}}\_2\right)^T\left(\mathbf{B}\_3\tilde{\mathbf{A}}\_3\right)}{||\mathbf{B}\_3\tilde{\mathbf{A}}\_3||}\right) \tag{10}$$

$$q\_{3,4} = \gamma \tag{11}$$

Here, **e**2 is the unit vector along **y**, and **~ R** denotes the initial state of **R**, which is given as follows:

$$\mathbf{R} = \begin{bmatrix} \cos \gamma & 0 & \sin \gamma \\ \sin \beta \sin \gamma & \cos \beta & -\sin \beta \cos \gamma \\ -\cos \beta \sin \gamma & \sin \beta & \cos \beta \cos \gamma \end{bmatrix} \tag{12}$$

#### *3.2. Actual Forward Kinematics*

The nominal inverse kinematics described above does not consider geometric errors. However, geometric errors exist in the 2UPR-RPU parallel manipulator. In this section, the actual forward kinematics of the limbs in the manipulator is derived in detail.

As shown in Figure 2, four local coordinate systems {**F***i*,*j*; **x***i*,*j*, **y***i*,*j*, **z***i*,*j*} are assigned to each limb to describe the geometric errors of the 2UPR-RPU parallel manipulator, where the initial pose of the moving platform is - *z*0 *β*0 *γ*0 T = - −0.2m 0 0T under the home configuration. The coordinate systems {**o**B; **x**, **y**, **z**} and {**o**A; **u**, **v**, **w**} are identical to those in Figure 1. For brevity, we use {**F***i*,*j*} instead of {**F***i*,*j*; **x***i*,*j*, **y***i*,*j*, **z***i*,*j*}. It is worth mentioning that this figure only shows **x***i*,*j* and **z***i*,*j* of the local coordinate systems, and **y***i*,*j* can be determined according to the right-hand rule, which will not be illustrated in detail here. The definitions of the local coordinate systems for the two UPR limbs and the RPU limb are listed in Tables 1–3.

**Figure 2.** 2UPR-RPU parallel manipulator and its local coordinate systems.

**Table 1.** Definitions of local coordinate systems for the first UPR limb.

**Table 2.** Definitions of local coordinate systems for the second UPR limb.

**Table 3.** Definitions of local coordinate systems for the RPU limb.

The end poses of the *i*th limb can be obtained from the local product of the exponential formula [25] as

$$\mathbf{g}\_{i}(\mathbf{q}\_{i}) = \mathbf{g}\_{i,0}e^{\hat{\boldsymbol{\xi}}\_{i,1}q\_{i,1}}\mathbf{g}\_{i,1}e^{\hat{\boldsymbol{\xi}}\_{i,2}q\_{i,2}}\mathbf{g}\_{i,2}e^{\hat{\boldsymbol{\xi}}\_{i,3}q\_{i,3}}\mathbf{g}\_{i,3}e^{\hat{\boldsymbol{\xi}}\_{i,4}q\_{i,4}}\mathbf{g}\_{i,4}, i = 1,2,3\tag{13}$$

where **g***i* denotes the homogeneous transformation matrix (HTM) of {**o**A; **u**, **v**, **w**} with respect to {**o**B; **x**, **y**, **z**} calculated using the *i*th limb. ζ*i,j* denotes the screw coordinates of **s***i,j* with respect to {**F***i*,*j*}, which can be written as [25,26]

$$\begin{cases} \mathcal{L}\_{i,j} = \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}^T \text{for R joint} \\ \mathcal{L}\_{i,j} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}^T \text{for P joint} \end{cases} \tag{14}$$

Here, *e* **^** ζ*i*,*jqi*,*j* denotes the exponential map from the Lie algebra *se*(3) to the special Euclidean group *SE*(3), which can be obtained using (A1)–(A4) in Appendix A. **g***i*,j is the HTM between adjacent coordinate systems when the parallel manipulator is under the home configuration. To be more specific, **g***i*,0 denotes the HTM of {**F***i*,1} with respect to {**o**B; **x**, **y**, **z**}; **g***i*,4 denotes the HTM of {**o**A; **u**, **v**, **w**} with respect to {**F***i*,4}; when *j* = 0 and *j* = 4, **g***i*,*j* is the HTM of {**F***i*,*j*+1} with respect to {**F***i*,*j*}. **g***i*,*j* can be written as

$$\begin{cases} \mathbf{g}\_{1,0} = \mathbf{Trans}(x, l\_{\mathrm{B}} + d) \mathbf{Rot}(y, \pi/2) \\ \mathbf{g}\_{1,1} = \mathbf{Trans}(z, -d) \mathbf{Rot}(x, -\pi/2) \mathbf{Rot}(z, \tilde{q}\_{1,2} - \pi/2) \\ \mathbf{g}\_{1,2} = \mathbf{Rot}(y, \pi/2) \\ \mathbf{g}\_{1,3} = \mathbf{Trans}(z, \tilde{q}\_{1,3}) \mathbf{Rot}(y, -\pi/2) \mathbf{Rot}(z, -\tilde{q}\_{1,2} + \pi/2) \\ \mathbf{g}\_{1,4} = \mathbf{Trans}(y, l\_{\mathrm{A}}) \mathbf{Rot}(y, -\pi/2) \mathbf{Rot}(z, -\pi/2) \end{cases} \tag{15}$$

$$\begin{cases} \mathbf{g}\_{2,0} = \mathsf{Trans}(x, -l\_{\mathrm{B}} - d)\mathbf{Rot}(y, \pi/2) \\ \mathbf{g}\_{2,1} = \mathsf{Trans}(z, d)\mathbf{Rot}(x, -\pi/2)\mathbf{Rot}(z, \tilde{q}\_{2,2} - \pi/2) \\ \mathbf{g}\_{2,2} = \mathsf{Rot}(y, \pi/2) \\ \mathbf{g}\_{2,3} = \mathsf{Trans}(z, \tilde{q}\_{2,3})\mathbf{Rot}(y, -\pi/2)\mathbf{Rot}(z, -\tilde{q}\_{2,2} + \pi/2) \\ \mathbf{g}\_{2,4} = \mathsf{Trans}(y, -l\_{\mathrm{A}})\mathbf{Rot}(y, -\pi/2)\mathbf{Rot}(z, -\pi/2) \end{cases} \tag{16}$$

$$\begin{cases} \mathbf{g}\_{3,0} = \mathsf{Trans}(y, l\_{\mathrm{B}})\mathbf{Rot}(y, \pi/2)\mathbf{Rot}(z, \tilde{q}\_{3,1} - \pi/2) \\ \mathbf{g}\_{3,1} = \mathsf{Rot}(y, \pi/2) \\ \mathbf{g}\_{3,2} = \mathsf{Trans}(z, \tilde{q}\_{3,2})\mathbf{Rot}(y, -\pi/2)\mathbf{Rot}(z, -\tilde{q}\_{3,1} + \pi/2) \\ \mathbf{g}\_{3,3} = \mathsf{Trans}(y, -c)\mathbf{Rot}(x, -\pi/2) \\ \mathbf{g}\_{3,4} = \mathsf{trans}(z, c - l\_{\mathrm{A}})\mathbf{Rot}(y, -\pi/2)\mathbf{Rot}(z, -\pi/2) \end{cases} \tag{17}$$

where **Trans**(*x*, *l*B) denotes the HTM that translates by *l*B along *x*, and **Rot**(*y*, *π*/2) denotes the HTM that rotates by *π*/2 around *y*. *q*1,2 is the initial angle between **x** and **B**1**A**1, *q*2,2 is the initial angle between **B**2**B**1 and **B**2**A**2, and *q*3,1 is the initial angle between **B**3**o**B and **B**3**A**3, which can be expressed as

$$\tilde{q}\_{1,2} = \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_1 \tilde{\mathbf{A}}\_1}{||\mathbf{B}\_1 \tilde{\mathbf{A}}\_1||}\right) \tag{18}$$

$$\tilde{q}\_{2,2} = \arccos\left(\frac{\mathbf{e}\_1^T \mathbf{B}\_2 \tilde{\mathbf{A}}\_2}{||\mathbf{B}\_2 \tilde{\mathbf{A}}\_2||}\right) \tag{19}$$

$$\tilde{q}\_{3,1} = \arccos\left(\frac{-\mathbf{e}\_2^T \mathbf{B}\_3 \tilde{\mathbf{A}}\_3}{||\mathbf{B}\_3 \tilde{\mathbf{A}}\_3||}\right) \tag{20}$$

In contrast to *q*1,2, *q*2,2, and *q*3,1, *q*1,3, *q*2,3, and *q*3,2 are the initial positions of the actuated P-joints, and we have

$$
\widetilde{q}\_{1,3} = \|\mathbf{B}\_1 \widetilde{\mathbf{A}}\_1\| \tag{21}
$$

$$\tilde{q}\_{2,3} = \|\mathbf{B}\_2 \tilde{\mathbf{A}}\_2\|\tag{22}$$

$$
\tilde{q}\_{3,2} = ||\mathbf{B}\_3 \tilde{\mathbf{A}}\_3|| \tag{23}
$$

The linear errors δ*i*,*j* of {**F***i*,*j*+1} along **x***i*,*j*, **y***i*,*j*, and **z***i*,*j* can be expressed as follows:

$$\mathfrak{G}\_{i,j} = \begin{bmatrix} \delta^x\_{i,j} & \delta^y\_{i,j} & \delta^z\_{i,j} \end{bmatrix}^\mathrm{T}, i = 1, 2, 3 \text{ and } j = 0, \cdots, 4 \tag{24}$$

In addition to linear errors, angular errors also exist. The angular errors ε*i*,*j* of {**F***i*,*j*+1} around **x***i*,*j*, **y***i*,*j*, and **z***i*,*j* can be expressed as follows:

$$\mathfrak{e}\_{i,j} = \begin{bmatrix} \mathfrak{e}\_{i,j}^{x} & \mathfrak{e}\_{i,j}^{y} & \mathfrak{e}\_{i,j}^{z} \end{bmatrix}^{\mathrm{T}}, i = 1, 2, 3 \text{ and } j = 0, \cdots, 4 \tag{25}$$

where δ*i*,0 and ε*i*,0 denote the linear and angular errors of {**F***i*,1} with respect to {**o**B; **x**, **y**, **z**}, respectively. δ*i*,4 and ε*i*,4 denote the linear and angular errors of {**F***i*,4} with respect to {**o**A; **u**, **v**, **w**}, respectively. Among the 90 error parameters, *εx* 1,0, *ε y* 1,1, *δz* 1,3, *εx* 1,3, *εx* 2,0, *ε y* 2,1, *δz* 2,3, *εx* 2,3, *εx* 3,0, *δz* 3,2, *εx* 3,2, and *ε y* 3,3 represent the initial displacement errors of the 12 joints. In addition, the values of *δx* 1,2, *δx* 2,2, *ε y* 1,4, *ε y* 2,4, *δx* 3,1, *εz* 3,3, and *ε y* 3,4 are zeros since the definitions of local coordinate systems. Therefore, the rest 71 error parameters represent the linear and angular geometric errors.

Setting the values of error parameters other than geometric errors to zeros, the HTM of the geometric errors between adjacent coordinate systems can be written as

$$
\Delta \mathbf{g}\_{i,j} = \begin{bmatrix}
\hat{\mathfrak{e}}\_{i,j} & \mathfrak{S}\_{i,j} \\
0\_{1\times3} & 1
\end{bmatrix}, i = 1, 2, 3 \text{ and } j = 0, \cdots, 4 \tag{26}
$$

where *e* **^** ε*i*,*j* denotes the exponential map from the Lie algebra *so*(3) to the special orthogonal group *SO*(3), which can be determined using (A3) and (A5) in Appendix A.

The end poses of the *i*th limb that include sthe linear and angular geometric errors can then be obtained as follows:

$$\mathbf{g}\_{i}^{\mathcal{J}\varepsilon}(\mathbf{q}\_{i}) = \Delta \mathbf{g}\_{i,0} \mathbf{g}\_{i,0} \epsilon^{\hat{\mathcal{L}}\_{i,1}\mathbf{q}\_{i,1}} \Delta \mathbf{g}\_{i,1} \mathbf{g}\_{i,1} \epsilon^{\hat{\mathcal{L}}\_{i,2}\mathbf{q}\_{i,2}} \Delta \mathbf{g}\_{i,2} \mathbf{g}\_{i,2} \epsilon^{\hat{\mathcal{L}}\_{i,3}\mathbf{q}\_{i,3}} \Delta \mathbf{g}\_{i,3} \mathbf{g}\_{i,3} \epsilon^{\hat{\mathcal{L}}\_{i,4}\mathbf{q}\_{i,4}} \mathbf{g}\_{i,4} \Delta \mathbf{g}\_{i,4}^{-1}, i = 1,2,3 \tag{27}$$

which can be rewritten as

$$\mathbf{g}\_{i}^{\circ \varepsilon}(\mathbf{q}\_{i}) = \mathbf{g}\_{i,0}^{\circ \varepsilon} e^{\hat{\mathcal{L}}\_{i,1}q\_{i,1}} \mathbf{g}\_{i,1}^{\varepsilon \varepsilon} e^{\hat{\mathcal{L}}\_{i,2}q\_{i,2}} \mathbf{g}\_{i,2}^{\varepsilon \varepsilon} e^{\hat{\mathcal{L}}\_{i,3}q\_{i,3}} \mathbf{g}\_{i,3}^{\varepsilon \varepsilon} e^{\hat{\mathcal{L}}\_{i,4}q\_{i,4}} \mathbf{g}\_{i,4}^{\varepsilon \varepsilon}, i = 1,2,3\tag{28}$$

#### **4. Evaluation Model of Deformations**

As mentioned previously, the 2UPR-RPU parallel manipulator is over-constrained. Theoretically, the end poses of any two limbs can also be consistent with each other through the motion deviations of passive joints when the internal-force-and-deformation-related geometric errors are zero, which can be expressed as

$$\mathbf{g}\_{i}^{\mathcal{K}} + \Delta \mathbf{g}\_{i}^{\mathcal{K}} = \mathbf{g}\_{k}^{\mathcal{K}} + \Delta \mathbf{g}\_{k}^{\mathcal{K}} \tag{29}$$

where Δ**g***ge i* and Δ**g***ge k* denote the end-pose deviations of the *i*th and *k*th limbs caused by the motion deviations of passive joints, respectively. The end-pose deviation between the *i*th and *k*th limbs can be written as [25,26]

$$
\Delta \mu\_{k,i} = \left\{ \log \left[ \mathbf{g}\_k^{\mathbb{S}^c} \left( \mathbf{g}\_i^{\mathbb{S}^c} \right)^{-1} \right] \right\}^\vee \tag{30}
$$

where log[·] stands for the logarithmic operation from *SE*(3) to *se*(3), and it can be obtained using (A6) and (A7) in Appendix A. ∨ represents the reverse operation of (A1). The end-pose deviation can be rewritten in screw form as follows:

$$
\Delta\mu\_{k,\bar{l}} = \Delta\mu\_{\bar{l}} - \Delta\mu\_k \tag{31}
$$

where the screws Δμ*i* and Δμ*k* denote the end-pose deviations of the *i*th and *k*th limbs originating from the motion deviations of passive joints, respectively. Take Δμ*i* as an example. Taking the partial differential of (28) with respect to the displacements of the passive joints, Δμ*i* can be expressed as follows:

$$
\Delta \mu\_i = \Psi\_i \Phi\_i \Delta \mathbf{q}\_{i\prime} i = 1,2,3 \tag{32}
$$

where Δ**q***i* denotes the motion deviation of the passive joints of the *i*th limb.

When *i* = 1 and *i* = 2, we have

$$
\Delta \mathbf{q}\_i = \begin{bmatrix} \Delta q\_{i,1} & \Delta q\_{i,2} & \Delta q\_{i,4} \end{bmatrix}^T \tag{33}
$$

and when *i* = 3, we have

$$
\Delta \mathbf{q}\_i = \begin{bmatrix} \Delta q\_{i,1} & \Delta q\_{i,3} & \Delta q\_{i,4} \end{bmatrix}^T \tag{34}
$$

For the coefficient matrices **Ψ***i* and **Φ***i*, when *i* = 1 and *i* = 2, we obtain

$$\mathbf{Y}\_i = \begin{bmatrix} \mathbf{I}\_6 & \operatorname{Ad}\left(e^{\hat{\tilde{\mathbf{k}}}\_{i,1}q\_{i,1}}\right) & \operatorname{Ad}\left(e^{\hat{\tilde{\mathbf{k}}}\_{i,1}q\_{i,1}}e^{\hat{\tilde{\mathbf{k}}}\_{i,2}q\_{i,2}}e^{\hat{\tilde{\mathbf{k}}}\_{i,3}q\_{i,3}}\right) \end{bmatrix} \in \mathbb{R}^{6 \times 18} \tag{35}$$

$$\Phi\_i = \text{Blkdiag}(\mathfrak{E}\_{i,1'}, \mathfrak{E}\_{i,2'}, \mathfrak{E}\_{i,4}) \in \mathbb{R}^{18 \times 3} \tag{36}$$

When *i* = 3, we obtain

$$\mathbf{V}\_{i} = \begin{bmatrix} \mathbf{I}\_{6} & \operatorname{Ad}\left(e^{\hat{\underline{\mathbf{k}}}\_{i,1}q\_{i,1}}e^{\hat{\underline{\mathbf{k}}}\_{i,2}q\_{i,2}}\right) & \operatorname{Ad}\left(e^{\hat{\underline{\mathbf{k}}}\_{i,1}q\_{i,1}}e^{\hat{\underline{\mathbf{k}}}\_{i,2}q\_{i,2}}e^{\hat{\underline{\mathbf{k}}}\_{i,3}q\_{i,3}}\right) \end{bmatrix} \in \mathbb{R}^{6 \times 18} \tag{37}$$

$$\Phi\_i = \text{Blkdiag}\left(\mathfrak{E}\_{i,1}, \mathfrak{E}\_{i,3}, \mathfrak{E}\_{i,4}\right) \in \mathbb{R}^{18 \times 3} \tag{38}$$

where **I**6 is an identity matrix of order six. Ad(·) is an adjoint representation of *SE*(3) and is given in (A8) in Appendix A. Blkdiag(·) denotes a block-diagonal matrix. ξ*i,j* denotes the screw coordinates of **s***i,j* with respect to {**o**B; **x**, **y**, **z**}, which can be written as follows [25]:

$$\mathfrak{E}\_{i,j} = \text{Ad}\left(\mathbf{g}\_{i,0'}^{\xi c} \mathbf{g}\_{i,1'}^{\xi c} \cdot \cdots \mathbf{g}\_{i,j-1}^{\xi c}\right) \mathfrak{L}\_{i,j} \tag{39}$$

Combining (31) with (32) yields

$$
\begin{bmatrix}
\Delta\mu\_{2,1} \\
\Delta\mu\_{3,2} \\
\Delta\mu\_{1,3}
\end{bmatrix} = \begin{bmatrix}
\Psi\_1\Phi\_1 & -\Psi\_2\Phi\_2 & 0 \\
0 & \Psi\_2\Phi\_2 & -\Psi\_3\Phi\_3 \\
\end{bmatrix} \begin{bmatrix}
\Delta\mathbf{q}\_1 \\
\Delta\mathbf{q}\_2 \\
\Delta\mathbf{q}\_3
\end{bmatrix} \tag{40}
$$

Let

$$
\Delta\boldsymbol{\mu} = \begin{bmatrix} \Delta\boldsymbol{\mu}\_{2,1}^{\mathrm{T}} & \Delta\boldsymbol{\mu}\_{3,2}^{\mathrm{T}} & \Delta\boldsymbol{\mu}\_{1,3}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}} \in \mathbb{R}^{18 \times 1} \tag{41}
$$

$$
\Delta \mathbf{q} = \begin{bmatrix} \Delta \mathbf{q}\_1^\mathrm{T} & \Delta \mathbf{q}\_2^\mathrm{T} & \Delta \mathbf{q}\_3^\mathrm{T} \end{bmatrix}^\mathrm{T} \in \mathbb{R}^{9 \times 1} \tag{42}
$$

$$\mathbf{J} = \begin{bmatrix} \Psi\_1 \Phi\_1 & -\Psi\_2 \Phi\_2 & 0 \\ 0 & \Psi\_2 \Phi\_2 & -\Psi\_3 \Phi\_3 \\ -\Psi\_1 \Phi\_1 & 0 & \Psi\_3 \Phi\_3 \end{bmatrix} \in \mathbb{R}^{18 \times 9} \tag{43}$$

Note that when Δμ is obtained using (30) and (41), the motion deviations of passive joints can be calculated as

$$
\Delta \mathbf{q} = \left(\mathbf{J}^{\mathrm{T}} \mathbf{J}\right)^{-1} \mathbf{J}^{\mathrm{T}} \Delta \boldsymbol{\mu} = \mathbf{J}\_c \Delta \boldsymbol{\mu} \tag{44}
$$

Based on the above work, an iterative model can be proposed to evaluate the deformations caused by geometric errors of the 2UPR-RPU over-constrained manipulator. The detailed processes are described below.

As shown in Figure 3, the proposed evaluation model mainly includes the following steps. Firstly, a target pose of the moving platform is input, and specified values are assigned to some of the 71 geometric errors; secondly, the nominal displacements of all joints are calculated based on the inverse kinematics; thirdly, the displacements of the passive joints are iteratively updated starting with the nominal values and the end condition is given as the maximum number of iterations or the target value of the infinity norm -Δμ*j* − Δμ*j*−1-∞; finally, the latest end-pose deviation Δμ*j* for the target pose is output. When the internal-force-and-deformation-related geometric errors are not all zeros, the end poses of the limbs cannot be consistent without deformations. Therefore, the latest Δμ*j* and indices based on it can be used to indirectly evaluate the limbs' comprehensive deformations caused by geometric errors of the 2UPR-RPU over-constrained manipulator.

**Figure 3.** Scheme I: Evaluation of the limbs' comprehensive deformations caused by geometric errors of the 2UPR-RPU over-constrained manipulator.

Considering that a large amount of matrix calculation is included in the proposed evaluation model, MATLAB is used for programming in Sections 5 and 6.

#### **5. Geometric Error Identification**

Finding the internal-force-and-deformation-related geometric errors is the basis of sensitivity analysis. In this section, the reachable workspace of the 2UPR-RPU parallel manipulator is described. For geometric error identification and verification, 692 and 1738 target poses are selected in the reachable workspace. Subsequently, internal-forceand-deformation-related geometric errors in the manipulator are identified based on the proposed evaluation model and an evaluation index. Finally, simulations are conducted to verify the correctness of the identification results.

#### *5.1. Identification Analysis*

The structural parameters of the 2UPR-RPU parallel manipulator are presented in Table 4. Using the space search method [29], the reachable workspace of the manipulator can be obtained. The search results are shown in Figure 4. Because the end poses at the boundaries of the reachable workspace are more sensitive to geometric errors, the 692 target poses shown in Figure 5 are uniformly selected for geometric error identification. To identify the internal-force-and-deformation-related geometric errors, the evaluation index of the maximum comprehensive deformation of a limb can be written as

$$
\Delta\mu\_{\text{max}} = \max\left( \|\Delta\mu\_{2,1}^j\|\_{\prime} \|\Delta\mu\_{3,2}^j\|\_{\prime} \|\Delta\mu\_{1,3}^j\| \right) \tag{45}
$$

**Table 4.** Structural parameters of the 2UPR-RPU parallel manipulator.

**Figure 4.** Reachable workspace of the 2UPR-RPU parallel manipulator.

**Figure 5.** 692 target poses of the 2UPR-RPU parallel manipulator.

Based on Scheme I and (45), 692 Δ*μ*maxs can be calculated for the selected target poses of the moving platform. If the 692 Δ*μ*maxs are not close to zero, it means that there are internal-force-and-deformation-related geometric errors among the geometric errors that were assigned specified values. Without loss of generality, three groups of specified values for geometric errors are given, as listed in Table 5. The maximum iteration number *λ* and specified tolerance *τ* in Scheme I are set to 50 and 10−15, respectively.

**Table 5.** Specified geometric errors [13,26] for geometric error identification.

Taking group 1 as an example, the detailed processes are described as follows: (1) *δx* 1,0 is set to 0.005 m, and the remaining geometric errors are set to 0. (2) 692 Δ*μ*maxs are calculated according to Scheme I and (45). (3) If the number of Δ*μ*maxs that are smaller than 10−15 is less than 657 (≈95% of 692), then *δx* 1,0 is referred to as an internalforce-and-deformation-related geometric error. After repeating the above steps for the 71 geometric errors, 39 internal-force-and-deformation-related geometric errors were initially identified and are listed in Table 6. In the table, "-" denotes the internal-forceand-deformation-related geometric error; "–" denotes the error parameter that is not a geometric error.

**Table 6.** Initially identified internal-force-and-deformation-related geometric errors.

Some geometric errors between any two adjacent coordinate systems in a limb may be linearly dependent. Therefore, it is necessary to analyse geometric errors simultaneously. Based on the results in Table 6, the set of the six error parameters, *δx i*,*j* , *δ y i*,*j* , *δz i*,*j* ,*εx i*,*j* ,*ε y i*,*j* ,*εz i*,*j* , are regarded as one unit. Take *δx* 1,1, *δ y* 1,1, *δz* 1,1,*εx* 1,1,*ε y* 1,1,*εz* 1,1 as an example. *δx* 1,1, *δ y* 1,1, *δz* 1,1,*εx* 1,1,*ε y* 1,1,*εz* 1,1 is set to [0.005 m, 0, 0.005 m, 0, 0, 0.005 rad], and the remaining units are set to [0,0,0,0,0,0]. Then, 692 Δ*μ*maxs are calculated according to Scheme I and (45). If the number of Δ*μ*maxs that are smaller than 10−15 is less than 657, the internal-force-and-deformation-related geometric errors are included in *δx* 1,1, *δz* 1,1, and *εz* 1,1. Then, *δx* 1,1, *δz* 1,1, and *εz* 1,1 are set to 0 in turn, and the remaining units are unchanged. The Δ*μ*maxs are recalculated. If the number of Δ*μ*maxs that decrease significantly is greater than 656, it is determined that the geometric error, which is set as 0, will cause internal forces and deformations. These steps were repeated for each error unit and the results are listed in Table 7. The identification results for groups 2 and 3 in Table 5 are the same as those shown in Table 7. The results demonstrate that there are 41 internal-force-anddeformation-related geometric errors, where the number of angular geometric errors is greater than that of linear geometric errors. In addition, the internal-force-and-deformation-related geometric errors of the first UPR limb are the same as those of the second UPR limb because of

the symmetric distribution of the two limbs. For the RPU limb, the geometric errors that cause internal forces and deformations are angular geometric errors.

**Table 7.** Identified internal-force-and-deformation-related geometric errors.

#### *5.2. Simulation Analysis*

To validate the correctness of the identified results listed in Table 7, three groups of numerical simulations were conducted using 1738 target poses of the 2UPR-RPU parallel manipulator, as shown in Figure 6. It is assumed that geometric errors are normally distributed with zero means [19,22]. Three groups of standard deviations are listed in Table 8. In the simulation, the internal-force-and-deformation-related geometric errors identified in Table 7 were set to 0, and the remaining 30 geometric errors were assigned random values generated by randn function using the standard deviations of δ*i*,*j* and ε*i*,*j* listed in Table 8. Then, according to Scheme I and (45), 1738 Δ*μ*maxs were calculated for each group. The simulation results are shown in Figure 7. It can be seen that Δ*μ*maxs of Group 1, Group 2, and Group 3, are all smaller than 10−15. This demonstrates that the internal-force-and-deformation-related geometric errors identified in Section 5.1 are correct.

**Figure 6.** 1738 target poses of the 2UPR-RPU parallel manipulator.

**Table 8.** Standard deviations of the geometric errors for the numerical simulations.

**Figure 7.** Simulation results using the standard deviations listed in Table 8. (**a**) Group 1; (**b**) Group 2; (**c**) Group 3.

#### **6. Sensitivity Analysis**

Sensitivity analysis can help reveal the influence of different internal-force-and-deformationrelated geometric errors on the limbs' comprehensive deformations. Since Δμ*j* is calculated iteratively in Scheme I, the Monte Carlo method [22] is utilised to conduct sensitivity analysis in this section. Two global sensitivity indices are proposed and the results of sensitivity analysis are verified through simulations.

#### *6.1. Sensitivity Indices*

According to (41), Δμ*j* consists of Δμ*j* 2,1, Δμ*j* 3,2, and Δμ*j* 1,3, which can be written as

$$
\Delta \boldsymbol{\mu}\_{k,i}^{j} = \begin{bmatrix} \boldsymbol{\omega}\_{k,i,1}^{j} & \boldsymbol{\omega}\_{k,i,2}^{j} & \boldsymbol{\omega}\_{k,i,3}^{j} & \boldsymbol{v}\_{k,i,1}^{j} & \boldsymbol{v}\_{k,i,2}^{j} & \boldsymbol{v}\_{k,i,3}^{j} \end{bmatrix}^{\mathrm{T}} \tag{46}$$

The end-orientation and end-position volumetric deviations between any two limbs are

$$
\Delta\omega\_{k,i}^{j} = \sqrt{\left(\omega\_{k,i,1}^{j}\right)^{2} + \left(\omega\_{k,i,2}^{j}\right)^{2} + \left(\omega\_{k,i,3}^{j}\right)^{2}}\tag{47}
$$

and

$$
\Delta v\_{k,i}^{j} = \sqrt{\left(v\_{k,i,1}^{j}\right)^2 + \left(v\_{k,i,2}^{j}\right)^2 + \left(v\_{k,i,3}^{j}\right)^2} \tag{48}
$$

Then, the evaluation indices of the average angular and linear comprehensive deformations of the three limbs can be written as

$$
\Delta\omega\_a^{\dot{j}} = \frac{\Delta\omega\_{2,1}^{\dot{j}} + \Delta\omega\_{3,2}^{\dot{j}} + \Delta\omega\_{1,3}^{\dot{j}}}{3} \tag{49}
$$

and

$$
\Delta v\_a^j = \frac{\Delta v\_{2,1}^j + \Delta v\_{3,2}^j + \Delta v\_{1,3}^j}{3} \tag{50}
$$

Under the condition that geometric errors are normally distributed with zero means, the sensitivity indices of the average angular and linear comprehensive deformations with respect to a geometric error can be written as

$$\mu\_{\omega,p} = \frac{\sigma\left(\Delta\omega\_{a,p}^j\right)}{\sigma\left(G\varepsilon\_p\right)}, \ p = 1, 2, \cdots, \ 25\tag{51}$$

and

$$\mu\_{v,p} = \frac{\sigma\left(\Delta v\_{a,p}^{j}\right)}{\sigma\left(G\varepsilon\_{p}\right)}, \ p = 1, 2, \cdots, 25\tag{52}$$

where *σ*(·) denotes the standard deviation and can be calculated by **std** function. *Gep*s are the internal-force-and-deformation-related geometric errors of the first and third limbs. Because of the symmetric distribution of the first and second limbs, the internal-force-anddeformation-related geometric errors of the second limb are not considered. Generally, the values of sensitivity indices vary with different target poses of the moving platform. Hence, *m* target poses should be chosen and the global sensitivity indices can be written as [16]

$$
\mu\_{\omega,p}^{\mathcal{S}} = \frac{\sum\_{i=1}^{m} \mu\_{\omega,p,i}}{m} + \sigma(\mu\_{\omega,p}) \tag{53}
$$

and

$$
\mu\_{\upsilon,p}^{\mathcal{S}} = \frac{\sum\_{i=1}^{m} \mu\_{\upsilon,p,i}}{m} + \sigma \left(\mu\_{\upsilon,p}\right) \tag{54}
$$

#### *6.2. Sensitivity Analysis*

Based on the equations in Section 6.1 and Scheme I, the detailed processes to calculate the two global sensitivity indices with respect to *Gep* are described as follows: (1) Set *σ* # *Gep* \$ to 1 mm (0.001 m) or 1◦ (π/180 rad) for linear or angular geometric error. And the other 40 internal-force-and-deformation-related geometric errors are set to 0. In addition, the rest 30 linear or angular geometric errors are set to 1 mm or 1◦. (2) Assign 1000 random values that obey the normal distribution to *Gep* and calculate 1000 Δ*ωj a*,*p*s and Δ*v j a*,*p*s. (3) Calculate *σ* Δ*ωj a*,*p* , *μω*,*p*, *σ* Δ*v j a*,*p* , and *μv*,*p* for a target pose of the moving platform. (4) Repeat the above steps for *m* target poses and calculate the global sensitivity indices*μg ω*,*p* and *μg v*,*p*.

In order to improve the computational efficiency, 158 of the 1738 target poses shown in Figure 6 were selected uniformly to perform the above steps for each *Gep*. The global sensitivity indices of the average angular and linear comprehensive deformations with respect to *Gep*s are shown in Figures 8 and 9, respectively. It can be seen that the values of *μg ω*,*p* with respect to *Ge*5(*δ y* 1,1), *Ge*7(*εz* 1,1), *Ge*8(*δz* 1,2), *Ge*11(*δx* 1,3), and *Ge*14(*δ y* 1,4), are zero. This indicates that the corresponding geometric errors have no effects on the average angular comprehensive deformation. It is worth mentioning that *δ y* 2,1, *εz* 2,1, *δz* 2,2, *δx* 2,3, and *δ y* 2,4, have also no effects on the average angular comprehensive deformation due to the symmetric distribution of the first and second limbs. Comparing Figure 8 with Figure 9, it can also be found that the value of *μg v*,*p* is larger than that of *μg ω*,*p* for each *Gep*. This demonstrates that the internal-force-and-deformation-related geometric errors have greater effects on the average linear comprehensive deformation. Thus, the distribution of the global sensitivity index *μg v*,*p* is more useful for accuracy synthesis. According to Figure 9, *Gep*s can be sorted in descending order as follows: *Ge*16(*εz* 1,4), *Ge*25(*εz* 3,4), *Ge*15(*εx* 1,4), *Ge*4(*εz* 1,0), *Ge*24(*εx* 3,4), *Ge*12(*ε y* 1,3), *Ge*6(*εx* 1,1), *Ge*3(*ε y* 1,0), *Ge*13(*εz* 1,3), *Ge*9(*εx* 1,2), *Ge*18(*εz* 3,0), *Ge*23(*εx* 3,3), *Ge*17(*ε y* 3,0), *Ge*19(*εx* 3,1), *Ge*22(*εz* 3,2), *Ge*10(*ε y* 1,2), *Ge*20(*ε y* 3,1), *Ge*21(*ε y* 3,2), *Ge*8(*δz* 1,2), *Ge*14(*δ y* 1,4), *Ge*5(*δ y* 1,1), *Ge*11(*δx* 1,3), *Ge*1(*δ y* 1,0), *Ge*2(*δz* 1,0), and *Ge*7(*εz* 1,1). In order to lower the cost of fabrication and assembly, the allowable range of geometric errors should be larger and larger from *Ge*16 to *Ge*7.

**Figure 8.** Global sensitivity of the average angular comprehensive deformation with respect to *Gep*s.

#### *6.3. Verification*

6.3.1. Average Angular Comprehensive Deformation

As shown in Table 9, three groups of specified values for geometric errors are given. For each group, 1738 Δ*ωj a*s were calculated according to Scheme I and using the target poses shown in Figure 6. The maximum and average values of Δ*ωj a* are listed in Table 9. It can be seen that both the maximum and average values of Δ*ωj a* do not change from Group 1 to Group 3. This indicates that *Ge*5, *Ge*7, *Ge*8, *Ge*11, *Ge*14, *δ y* 2,1, *εz* 2,1, *δz* 2,2, *δx* 2,3, and *δ y* 2,4, have no effects on the average angular comprehensive deformation.

**Table 9.** Sensitivity analysis results of the average angular comprehensive deformation.

#### 6.3.2. Average Linear Comprehensive Deformation

As shown in Table 10, *Ge*7(*εz* 1,1), which has the smallest effect on the average linear comprehensive deformation, is set to 1◦. Then, the other *Gep*s are set to *μg v*,7/*μg v*,*p* mm or ◦. It is worth mentioning that the corresponding internal-force-and-deformation-related geometric errors of the second limb are assigned the same values as the first limb due to the symmetric distribution of the two UPR limbs. The remaining 30 geometric errors are set to 0.1 mm or ◦. According to Scheme I and using the target poses shown in Figure 6, 1738Δ*ωj a*s and Δ*v j a*s were calculated. For comparison, the internal-force-and-deformation-related linear and angular geometric errors are set to their average values, 0.0209 mm and 0.0743 ◦, respectively, while the values of the remaining 30 geometric errors are unchanged. After recalculation, the maximum and average values of Δ*ωj a* and Δ*v j a* are listed in Table 11. It can be seen that both the maximum and average values of Δ*v j a* are larger than that of Δ*ωj a* for each group. This indicates that the internal-force-and-deformation-related geometric errors have greater effects on the average linear comprehensive deformation. It can also be found that from Group 2 to Group 1, the maximum and average values of Δ*ωj a* and Δ*v j a* decreased by 84%, 83%, 91%, and 89%, respectively. This demonstrates that at the same cost, restricting the values of geometric errors according to the sensitivity analysis results of the average linear comprehensive deformation can dramatically decrease the average angular and linear comprehensive deformations. Furthermore, it indirectly verifies the sensitivity analysis results of the average linear comprehensive deformation.

**Table 10.** Specified geometric errors for verification.

**Table 11.** Sensitivity analysis results of the average linear comprehensive deformation.

#### **7. Conclusions**

This paper deals with error modelling and sensitivity analysis of geometric errors that cause internal forces and deformations in the 2UPR-RPU over-constrained parallel manipulator. Conclusions are drawn as follows:

(1) The nominal inverse kinematics and actual forward kinematics of the over-constrained parallel manipulator are analysed according to the vector theory and the local product of the exponential formula. On this basis, an iterative model is established to indirectly evaluate the limbs' comprehensive deformations caused by geometric errors.

(2) Based on the iterative evaluation model, the maximum Euclidean norm of the endpose deviations of limbs is defined as an evaluation index of the maximum comprehensive deformation of a limb. Programming with MATLAB, 41 internal-force-and-deformationrelated geometric errors are identified. Among the 41 geometric errors, the number of angular geometric errors is greater than that of linear geometric errors; the geometric errors of the first UPR limb are the same as those of the second UPR limb; the geometric errors of the RPU limb are all angular geometric errors. The correctness of the identification results is verified through simulations under the condition that geometric errors are normally distributed with zero means.

(3) The global sensitivity indices of the average angular and linear comprehensive deformations with respect to internal-force-and-deformation-related geometric errors are proposed and calculated based on the Monte Carlo method. The results of sensitivity analysis demonstrate that *δ y* 1,1, *εz* 1,1, *δz* 1,2, *δx* 1,3, *δ y* 1,4, *δ y* 2,1, *εz* 2,1, *δz* 2,2, *δx* 2,3, and *δ y* 2,4, have no effects on the average angular comprehensive deformation. Furthermore, the internalforce-and-deformation-related geometric errors have greater effects on the average linear comprehensive deformation. Therefore, the distribution of the global sensitivity index of the average linear comprehensive deformation with respect to geometric errors is more meaningful for accuracy synthesis. Finally, the results of sensitivity analysis are verified through simulations.

Based on the work presented in this paper, we will establish a model for accuracy synthesis and determine the tolerances of the fabrication and assembly of the manipulator in the future.

**Author Contributions:** Conceptualization, X.D. and J.Z.; methodology, X.D.; software, B.W.; validation, X.D., B.W. and J.Z.; formal analysis, B.W.; investigation, X.D.; resources, X.D.; data curation, B.W.; writing—original draft preparation, X.D.; writing—review and editing, J.Z.; visualization, X.D.; supervision, J.Z.; project administration, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by [National Natural Science Foundation of China] grant number [52275469] and [Open Fund of State Key Laboratory of Robotics and System] grant number [SKLRS-2021-KF-08].

**Data Availability Statement:** Not applicable.

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

#### **Appendix A. Lie Groups and Lie Algebras**

Some equations about Lie groups and Lie algebras [25,30,31] are introduced here so that this work can be clearly understood. For a screw ζ = - ωT νTT , the ∧ operation denotes

$$
\hat{\mathbf{U}} = \begin{bmatrix} \hat{\mathbf{u}} & \mathbf{v} \\ \mathbf{0}\_{1 \times 3} & \mathbf{0} \end{bmatrix} \in \text{sec}(\mathbf{3}) \tag{A1}
$$

The exponential map from the Lie algebra *se*(3) to the special Euclidean group *SE*(3) can be determined by

$$e^{\hat{\mathcal{L}}q} = \begin{bmatrix} \mathbf{I}\_3 + \sin q \hat{\boldsymbol{\omega}} + (1 - \cos q) \hat{\boldsymbol{\omega}}^2 & \mathbf{V} \mathbf{v} \\ \mathbf{0}\_{1 \times 3} & 1 \end{bmatrix} \in SE(3) \tag{A2}$$

where **I**3 is an identity matrix of order three. **^** ω and **V** are expressed as

$$
\stackrel{\circ}{\omega} = \begin{bmatrix} 0 & -\omega\_3 & \omega\_2 \\ \omega\_3 & 0 & -\omega\_1 \\ -\omega\_2 & \omega\_1 & 0 \end{bmatrix} \in \text{so}(3) \tag{A3}
$$

$$\mathbf{V} = q\mathbf{I}\mathbf{\hat{s}} + (1 - \cos q)\mathbf{\hat{w}} + (q - \sin q)\mathbf{\hat{w}}^2 \tag{A4}$$

The exponential map from the Lie algebra *so*(3) to the special orthogonal group *SO*(3) can be determined by

$$\boldsymbol{\epsilon}^{\hat{\mathbf{w}}} = \mathbf{I}\_3 + \frac{\sin \|\boldsymbol{\omega}\|}{\left\|\boldsymbol{\omega}\right\|} \hat{\boldsymbol{\omega}} + \frac{1 - \cos \|\boldsymbol{\omega}\|}{\left\|\boldsymbol{\omega}\right\|^2} \hat{\boldsymbol{\omega}}^2 \tag{A5}$$

For a HTM **g** ∈ *SE*(3), the Lie algebra *se*(3) can be obtained as

$$\log(\mathbf{g}) = \frac{1}{8} \csc^3 \frac{\theta}{2} \sec \frac{\theta}{2} \begin{bmatrix} \theta \cos 2\theta - \sin \theta \\ -\theta \cos \theta - 2\theta \cos 2\theta + \sin \theta + \sin 2\theta \\ 2\theta \cos \theta + \theta \cos 2\theta - \sin \theta - \sin 2\theta \\ -\theta \cos \theta + \sin \theta \end{bmatrix}^{\mathrm{T}} \begin{bmatrix} \mathbf{I}\_4 \\ \mathbf{g} \\ \mathbf{g}^2 \\ \mathbf{g}^3 \end{bmatrix} \tag{A6}$$

where

$$\theta = \arccos\left(\frac{\text{Tr}(\mathbf{g}) - 2}{2}\right), \theta \in (-\pi, \pi) \tag{A7}$$

The adjoint representation of **g** can be written as

$$\operatorname{Ad}(\mathbf{g}) = \operatorname{Ad}\left(\begin{bmatrix} \mathbf{R} & \mathbf{t} \\ 0\_{1\times3} & 1 \end{bmatrix}\right) = \begin{bmatrix} \mathbf{R} & 0\_{3\times3} \\ \mathbf{t}\mathbf{R} & \mathbf{R} \end{bmatrix} \in \mathbb{R}^{6\times6} \tag{A8}$$

#### **References**

## *Article* **Consistent Solution Strategy for Static Equilibrium Workspace and Trajectory Planning of Under-Constrained Cable-Driven Parallel and Planar Hybrid Robots**

**Qingjuan Duan \*, Quanli Zhao and Tianle Wang**

School of Mechano-Electronic Engineering, Xidian University, Xi'an 710071, China **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** This paper presents a consistent solution strategy for static equilibrium workspaces of different types of under-constrained robots. Considering the constraint conditions of cable force and taking the least squares error of the static equilibrium equation as the objective, the convex optimization solution is carried out, and the static equilibrium working space of the under-constrained system is obtained. A consistent solution strategy is applied to solve the static equilibrium workspaces of the cable-driven parallel and planar hybrid robots. The dynamic models are presented and introducing parameters that are applied to make the system stable for point-to-point movements. Based on this model, the traditional polynomial-based point-to-point trajectory planning algorithm is improved by adding unconstrained parameters to the kinematic law function. The constraints of the dynamics model are incorporated into the trajectory planning process to achieve point-to-point trajectory planning for the under-constrained cable-driven robots. Finally, under-constrained cabledriven parallel robots with three cables and planar hybrid robot with two cables are taken as examples to carry out numerical simulation. The final results show that the point-to-point trajectory planning algorithm introducing parameters is effective and feasible and can provide theoretical guidance for the design of subsequent under-constrained robots.

**Keywords:** static equilibrium workspace; under-constrained cable-driven robot; consistent solution strategy; trajectory planning

#### **1. Introduction**

The cable-driven robot is a mechanism that employs cables in place of rigid-body to control the end-effector pose.

The classification of cable-driven parallel robots (CDPRs) was introduced by Ming and Higuchi [1]. A cable-driven robot with *n* degrees-of-freedom and *m* cables can be classified into the following four categories: (1) *n* + 1 < *m*: These robots are referred to as redundantly restrained positioning mechanisms (RRPM), the static forces of the robot are generally undefined. (2) *n* + 1 = *m*: These robots are called completely restrained positioning mechanisms (CRPM). All degrees of freedom can be controlled through cables. (3) *m* = *n*: This type of robot is called incompletely restrained positioning mechanism (IRPM). When external forces such as gravity applied, the robot is fully constrained. It can withstand a limited range of wrenches. (4) *m* < *n*: The robot is under-constrained positioning mechanism (URPM) and in general cannot withstand arbitrary external forces and torques. Due to the under-constrained nature, these robots have one feasible solution for cable tensions and mostly works under gravity conditions. This classification method is also applicable to cable-driven planar hybrid robots.

According to the structure of cable-driven robot, it is generally divided into series mechanism and parallel mechanism.

Cable-driven robot has the characteristics of simple structure, flexibility, large workspace, low inertia, high load rate, etc. It has a wide range of applications, such as: Five-hundred-meter

**Citation:** Duan, Q.; Zhao, Q.; Wang, T. Consistent Solution Strategy for Static Equilibrium Workspace and Trajectory Planning of Under-Constrained Cable-Driven Parallel and Planar Hybrid Robots. *Machines* **2022**, *10*, 920. https:// doi.org/10.3390/machines10100920

Academic Editors: Dan Zhang, Zhufeng Shao and Stéphane Caro

Received: 18 July 2022 Accepted: 29 September 2022 Published: 10 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Aperture Spherical radio Telescope (FAST), wind tunnel test, rehabilitation training, sports photography, etc. [2–5], it has become a hot spot in robotics research recent years.

The under-constrained cable-driven robots (UCR), with few drivers and low cost, has its special purpose, attracting more and more scholars' research interest. Carricato et al. [6–8] from Italy studied the cable-driven parallel robots with less than six cables, provided a geometrico-static model, and assessed the stability of static equilibrium within the framework of a constrained optimization problem. Several examples are provided, concerning robots with a number of cables that range from 2 to 4. Berti et al. [9] proposed a method based on interval analysis to solve the positive geometric statics problem of an under-constrained cable-driven parallel robot, and find all possible equilibrium poses of the end effector under a given cable length. Liu Xin et al. [10] proposed a consistent algorithm for solving the workspace of a cable-driven parallel robot under different constraints. Fu Ying et al. [11] conducted a dynamic analysis on the cable-driven system with four cables and six degrees of freedom. Zhao Zhigang et al. [12] proposed a comprehensive algorithm by combining the least squares method and the Monte Carlo algorithm to solve the statically balanced workspace for the cable-driven system with multi-robots. Peng Y et al. [13] analyzed the reachable workspace for spatial 3-cable under-constrained suspended cable driven parallel robots. The above-mentioned literature put forwards the solution method of static equilibrium workspace for under-constrained parallel robots. Based on this, a consistent solution strategy for the static equilibrium workspace of both cable-driven parallel & planar hybrid robots is put forward in this paper.

Ida E et al. [14] proposed a rest-to-rest trajectory planning for underactuated cabledriven parallel robots. Barbazza L et al. [15] design and optimally control an underactuated cable-driven micro–macro robot. Shang Weiwei et al. [16] proposed a geometrical approach to plan trajectories that extend beyond the static equilibrium workspace (SEW) of the mechanism. Zi B et al. [17] studied an algebraic method-based point to point trajectory planning of an under constrained cable suspended parallel robot with variable angle and height cable mast. Shi P et al. [18] studied Dimensional synthesis of a gait rehabilitation cable-suspended robot on minimum 2-norm tensions. The above-mentioned literature put forwards the planning trajectories for under-constrained parallel robots. Based on this, a consistent solution strategies of point-to-point trajectory planning introducing parameters for both cable-driven parallel & planar hybrid robots is put forward in this paper.

Firstly, the static equilibrium equations of the under constrained parallel and planar hybrid robots are established. Then, the mathematical description of the static equilibrium workspace of the under-constrained cable-driven robots (UCR) that meets the constraints of the driving motor power and cable strength is given. The characteristics of the static equilibrium equation and dynamics equation are analyzed, and the consistent solution strategy for the static equilibrium workspace and point-to-point trajectory planning introducing parameters of different types of under-constrained robots are given.

The paper is organized as follows: Section 2 presents the model for under-constrained parallel robots. Section 3 provides the model of under-constraint cable-driven planar hybrid robot. Section 4 describes the consistent solution strategy for static equilibrium workspace. Section 5 illustrates some results of static equilibrium workspace. Section 6 puts forward the consistent solution strategies for point-to-point trajectory planning. Section 7 illustrates the simulation results of trajectory planning. Finally, Section 8 draws conclusions.

#### **2. Model of Under-Constrained Cable-Driven Parallel Robots**

The schematic structure of a cable-driven robot with *n* degrees-of-freedom and *m* cables is shown in Figure 1. Here *m<n*, it's an under-constrained CDPR. The moving platform is connected to the base through m cables, the *ith* cable (*i = 1,2,* ... *,m*) exits from the fixed base at point *Ai*, connected to the moving platform at point *Bi*. the cable length is *Li*. *Oxyz* is a Cartesian coordinate which fixed to the base, and *O x y z* is the Cartesian coordinate fixed to the moving platform.

**Figure 1.** The schematic structure of a cable-driven parallel robot (If *m<n*, it is an underconstrained CDPR.).

*OPO* = [*OxO OyO OzO* ] T is the centroid of end effector in the *Oxyz* frame. *O Bi* = [*O xBi O yBi O zBi*] T (*i* = 1, 2, ··· , *m*) is the vector connecting point *O'* to the point *Bi* in the *O x y z* frame. *O Ai* = [*OxAi OyAi OzAi*] T (*i* = 1, 2, ··· , *m*) is the fixed base at point *Ai*, in the *Oxyz* frame.

*ORO* represents the rotational matrix from frame *O x y z* to frame *Oxyz*, in which *α*, *β*, *γ* are *x*-*y*-*z* the Euler angles.

$$\begin{aligned}{}^{O}R\_{\alpha\gamma} &= \text{rot}(x,\alpha)\text{rot}(y,\beta)\text{rot}(z,\gamma) = \\ &\begin{bmatrix} \text{c\beta c\gamma} & -c\beta s\gamma & s\beta\\ \cos\gamma + s\alpha s\beta c\gamma & c\alpha c\gamma - s\alpha s\beta s\gamma & -c\beta s\alpha\\ s\alpha s\gamma - c\alpha s\beta c\gamma & s\alpha c\gamma + c\alpha s\beta s\gamma & c\alpha c\beta \end{bmatrix} \end{aligned} \tag{1}$$

where *c* represents cos, *s* represents sin.

*OBi* is the vector in the *Oxyz* frame.

$${}^{O}B\_{i} = {}^{O}R\_{O} \, {}^{O}B\_{i} + {}^{O}P\_{O} \, {}^{} \tag{2}$$

*OLi* is the vector connecting point *Bi* to point *Ai* in the *Oxyz* frame. *ei* is the unit vector of *OLi*.

$$x\_i = \frac{^O L\_i}{||^O L\_i||\_2} = \frac{^O A\_i - ^O B\_i}{||^O A\_i - ^O B\_i||\_2} \tag{3}$$

Thus, the dynamics equations for end effector are shown as follows:

$$JT + G = \begin{bmatrix} m\ddot{\mathbf{x}} \\ I\dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times I\boldsymbol{\omega} \end{bmatrix} \tag{4}$$

where *T* = [*t*1 *t*2 ··· *tm*] T, *ti*(*i* = 1, 2, ··· , *m*) is cable tensions act on the end effector. *J* = - ˆ*J*1 ˆ*J*2 ··· ˆ*Jm* is the construction matrix, ˆ*Ji* = - *ei OBi* × *ei* T , and *G* = [0 0 − *mg* 0 0] T . .. *x* = [*O* .. *xO O* .. *yO O*.. *zO* ] T is the acceleration of the end-effector centroid *O* in the world coordinate system, and *I* = *ORO IO ORO* , *IO* is the inertia tensor of the end-effector in the local coordinate system.

*ω* and . *ω* are the angular velocity and angular acceleration of the end-effector. *ω*, . *ω* and *ε*, the Euler angle, satisfy the following relationship:

$$
\omega = H(\varepsilon)\dot{\varepsilon} = \begin{bmatrix} 1 & 0 & s\beta \\ 0 & ca & -sac\beta \\ 0 & sa & cac\beta \end{bmatrix} \begin{bmatrix} \dot{\alpha} \\ \dot{\beta} \\ \dot{\gamma} \end{bmatrix} \tag{5}
$$

. *ω* = . *H*. *ε* + *H*.. *ε* (6)

Substituting Equations (5) and (6) into Equation (4) yields a general expression for the system dynamics equation [19], i.e.,

$$\begin{aligned} M(q)\ddot{q} - s(q, \dot{q}) - J(q)T &= 0\\ M(q) = \begin{bmatrix} mE\_3 & 0\\ 0 & I H \end{bmatrix} \rho s(q, \dot{q}) &= \mathbf{G} - \begin{bmatrix} 0\\ I \dot{H} \dot{\varepsilon} + \omega \times I \omega \end{bmatrix} \end{aligned} \tag{7}$$

In Equation (7), .. *q* = -.. *x* .. *ε* T is the end-effector acceleration. *M* is the mass matrix of the end-effector. *s* is the force vector for a collection of Coriolis-type forces, gravitational forces, and external loads, etc. *J*(*q*)*T* is the cable-tension to which the end-effector is subjected.

Additionally, the static equilibrium equation of the system can be expressed by

$$JT + G = 0\tag{8}$$

#### **3. Model of Under-Constrained Cable-Driven Planar Hybrid Robot**

The schematic structure of a cable-driven planar hybrid robot with planar *n* -link serial robot with *n* degrees-of-freedom and *m* cables is shown in Figure 2. All links are connected by revolute joints to form a planar multi-link mechanism. This definition can also be generalized to space model. Here, *m<n*, it's an under-constrained cable-driven planar hybrid robot. {0} is the base frame, {*i*} is the link frame *xiyi*, *i*∈{1, 2, ··· , *n*}. *θi*(*i* = 1, 2, ··· , *n*) is the angle between link frame {*i* − 1} and {*i*}. - *θ*1 *θ*2 ··· *θn* T is the Joint coordinates of the system.

0*Ci* is the position vector of the *ith* link centroid in the fixed frame {0}. It can be expressed by

$$\prescript{0}{}{\mathbf{C}}\_{i} = \begin{bmatrix} ^0 \mathbf{x}\_{\mathbf{C}\_i} \, ^0 \mathbf{y}\_{\mathbf{C}\_i} \end{bmatrix}^T = \begin{pmatrix} ^0 \mathbf{R}\_1 ^1 \mathbf{R}\_2 \cdots \stackrel{i-1}{\cdot} \mathbf{R}\_i \end{pmatrix} ^i \mathbf{C}\_i \tag{9}$$

where 0*R*1 = ⎡ ⎣ c*θ*1 −*sθ*1 0 *sθ*1 *cθ*1 0 0 01 ⎤ ⎦, *i Ri*+1 = ⎡ ⎣ c*θi*+1 −*sθi*+1 0 *sθi*+1 *cθi*+1 *li* 0 01 ⎤ ⎦, *i Ri*+1 represents the rota-

tional matrix from frame *i* + 1 to frame *i*, *i* = 1, 2, ··· , *n* − 1. 0*Bi j* is the *jth* cable position vector of the *ith* link in the fixed frame. It can be expressed by

$${}^{0}B\_{j}^{i} = \left[ {}^{0}\mathbf{x}\_{B\_{j}^{i}} {}^{0}y\_{B\_{j}^{i}} \right]^{\mathbf{T}} = \left( {}^{0}\mathbf{R}\_{1} {}^{1}\mathbf{R}\_{2} \cdot \cdots {}^{i-1}\mathbf{R}\_{i} \right) {}^{i}B\_{j}^{i} \tag{10}$$

where *i Bi j* = [*i xBi j i yBi j* ] T is the cable position vector of the *ith* link in the link frame {*i*}.

**Figure 2.** The schematic structure of a cable−driven planar hybrid robot with planar serial *n*−link with *n* degrees-of-freedom and *m* cables (If *m* < *n*, it's an under-constrained cable-driven planar hybrid robot.).

The generalized force is solved by Lagrange equation,

$$\begin{aligned} \tau &= \frac{d}{d\_t} \left( \frac{\partial L}{\partial \dot{\theta}} \right) - \frac{\partial L}{\partial \theta} \\\\ L &= E\_k - E\_p \end{aligned} \tag{11}$$

$$\tau = \frac{\partial \left(\frac{\partial \dot{\theta}}{\partial \dot{\theta}}\right)}{\partial \dot{\theta}} \ddot{\theta} + \frac{\partial \left(\frac{\partial \dot{\theta}}{\partial \dot{\theta}}\right)}{\partial \theta} \dot{\theta} - \frac{\partial L}{\partial \theta} \tag{12}$$

where *τ* is the generalized force, *L* is the Lagrange function. *θ* = - *θ*1 *θ*2 ··· *θn* T is the generalized coordinate. *Ek* and *Ep* are the kinetic energy and potential energy of the system, respectively.

Write Equation (12) as a generic expression,

$$\begin{aligned} M(\theta)\ddot{\theta} - s(\theta, \dot{\theta}) &= \tau \\ M(\theta) = \frac{\partial \left(\frac{\partial I}{\partial \theta}\right)}{\partial \theta}, s(\theta, \dot{\theta}) &= -\frac{\partial \left(\frac{\partial I}{\partial \theta}\right)}{\partial \theta} \dot{\theta} + \frac{\partial I}{\partial \theta} \end{aligned} \tag{13}$$

where *M*(*θ*) is the inertia matrix, and *s*(*θ*, . *θ*) is the gravity, centrifugal force, Coriolis force, etc.

The expression (14) can be obtained according to the principle of virtual work

$$\tau^{\mathrm{T}}\delta\theta = \sum\_{i=1}^{m} T\_{\mathrm{i}} \mathbf{e}\_{\mathrm{i}}^{\mathrm{T}} \delta(^{0}B\_{\mathrm{i}}) = \sum\_{i=1}^{m} T\_{\mathrm{i}} \mathbf{e}\_{\mathrm{i}}^{\mathrm{T}} \left(\frac{\partial^{0} B\_{\mathrm{i}}}{\partial\theta} \delta\theta\right) \tag{14}$$

where *m* represents the number of cables, *T* is the cable tension, *e* is the unit vector of the cable tension, and 0*Bi* is the position of the *ith* pulling point in the coordinate system.

From Equation (14), the relationship between generalized force and cable tension can be obtained

$$\boldsymbol{\pi} = \boldsymbol{J\_{\overline{v}}}^{\mathrm{T}} \boldsymbol{T} \tag{15}$$

where

*Jv* T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *e*1 T *∂*0*B*1 *∂θ*1 *e*1 T *∂*0*B*1 *∂θ*2 ··· *e*1 T *∂*0*B*1 *∂θn e*2 T *∂*0*B*2 *∂θ*1 *e*2 T *∂*0*B*2 *∂θ*2 ··· *e*2 T *∂*0*B*2 *∂θn* . . . . . . ... . . . *em*T *∂*0*Bm ∂θ*1 *em*T *∂*0*Bm ∂θ*2 ··· *em*T *∂*0*Bm ∂θn* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ T

where *Jv* is the pseudo-Jacobian of the constraint equations, *Jv* T is the *n* × *m* matrix, *T* is the *m* × 1 cable tension matrix.

Combining Equations (13) and (15), the system dynamics equation is obtained as follows. .. .

$$M(\theta)\theta - s(\theta, \theta) - J\_v^{\text{T}}(\theta)T = 0\tag{16}$$

Since the system is in static equilibrium, and the generalized velocity . *θ* = 0, Formula (11) can be simplified as

$$
\pi = \frac{\partial E\_p}{\partial \theta} \tag{17}
$$

where *Ep* = *n* ∑ *i*=1 (−*nig*)0*yCi* . Taking the origin of the coordinate system as the reference point of potential energy.

Combine Equations (15) and (17), yields

$$J\_{\mathbb{P}}{}^{\mathrm{T}}T = \frac{\partial E\_{\mathbb{P}}}{\partial \theta} \tag{18}$$

#### **4. Consistent Solution Strategy for Static Equilibrium Workspace**

Combine Equations (8) and (18), yields

$$J(X)T = \mathcal{W}(X)\tag{19}$$

where *W* is the external force.

#### *4.1. The Definition of Static Equilibrium Workspace*

Due to the limitation of motor torque and cable strength, the tension of the cable must be within a certain range. Therefore, for a robot with *n* degrees of freedom pulled by *m* cables, the mathematical description of its static equilibrium workspace is as follows:

$$X = \begin{bmatrix} x\_1 & x\_2 & \cdots & x\_n \end{bmatrix}^T$$

$$\exists t\_{\min} \le t\_i \le t\_{\max} (t\_{\min} > 0 \cap t\_{\max} > 0 \cap t\_{\max} > t\_{\min}, i = 1, 2, \cdots, m) \tag{20}$$

$$J(X)T = \mathcal{W}(X)$$

where *X* represents the generalized coordinate. *t*min, *t*max are the minimum and the maximum allowable tension of the cable, respectively. *T* = [*t*1 *t*2 ··· *tm*] T is the cable tension, and *X* belongs to the Static equilibrium workspace.

#### *4.2. Consistent Solution Strategy*

For a cable-driven robot with *n* degrees-of-freedom and *m* cables, here *m<n*, it's an Under-constrained Cable-driven Robot (UCR). The Static equilibrium workspace of the under-constrained cable-driven system is analyzed as follows:

In static equilibrium Equation (19), *J* is *n* × *m* matrix, *T* is *m* × 1 matrix. For the under-constrained cable-driven system, the solution of the equation is in the form of the least squares solution, and it can be expressed by *T* = *J*†*W*. When this solution is within the limit of cable force, it is considered that the pose (generalized coordinate) is the static equilibrium point satisfying Equation (20).

Then, the following inequality holds

$$\begin{array}{l} \|f(X)T - \mathcal{W}(X)\|\_{2} < \sigma \ast\\ T = f^{\dagger}(X)\mathcal{W}(X), t\_{\text{min}} \le t\_i \le t\_{\text{max}} \end{array} \tag{21}$$

where *σ*∗ is the error of the least square solution.

To sum up, the consistent solution strategy for the Static equilibrium workspace of the UCRs is as follows:

#### **5. Simulation Results of Static Equilibrium Workspace**

#### *5.1. 3*−*6 Under-Constrained Parallel Robot*

Solve the Static equilibrium workspace of 3−6 (3 cables pulling 6 degrees of freedom) under-constrained parallel robot. The mass of end effector is *m* = 0.935 kg. See Table 1 for specific parameters.

**Table 1.** 3−6 Configuration of cable-driven parallel robot.

*t*min = 0.1 N is the lower limit of cable tension, and *t*max = 10 N the upper limit of cable tension. *σ*∗ = 1 × 10−8 is the error. The controllable degree of freedom is *Xa* = - *OxO* , *O yO* , *O zO* , and the uncontrollable degree of freedom is *Xb* = [*α*, *β*, *γ*], the rotation angle of *x*-axis, *y*-axis, and *z*-axis, respectively. Considering the volume of end effector and cable, set the search interval is

$$\mathbf{S}\_{\mathbf{x}} \in (-0.6, 0.6), \mathbf{S}\_{\mathbf{y}} \in (0.3, 1.6), \mathbf{S}\_{\mathbf{z}} \in (0.27, 1.52)$$

The search step of *δx*, *δy*, *δz* is 0.05m, and the pose set **Q** to be searched is generated.

Matlab is used to solve the static equilibrium workspace of 3−6 cable-driven parallel robot.

It can be seen from Figure 3 the 3D view, *xoy* planar view and *xoz* planar view of the Static equilibrium workspace that the Static equilibrium workspace of the 3−6 underconstrained parallel robot is approximately a triangular prism, whose cross section is symmetrical along the *x*-axis, and there are two narrow workspace cracks in the lower half (*z* < 0.75 m) of the triangular prism. Compared with the system schematic diagram 1, this situation occurs because when the center of mass of the end actuator is low, the tension force of the three cables cannot be balanced with gravity, so the pose of the end actuator cannot be kept stable. The anchor points *O A*1 and *O A*3 are symmetrical along the *y*-axis, and the anchor point *O A*2 is on the axis, so that the cross section of its Static equilibrium workspace is symmetrical about the *y*-axis.

#### *5.2. 0*−*0*−*2 Under-Constrained Planar Hybrid Robot*

Solve the Static equilibrium workspace of the 0−0−2 under-constrained planar hybrid robot (no cable attached on the first and second links, and two cables attach on the third link). The Tables 2 and 3 show the specific parameters of 0−0−2 cable-driven planar hybrid robot system.

**Table 2.** 0−0−2 configuration of cable-driven planar hybrid robot.

**Table 3.** Link parameters.

The lower limit of the cable tension is *t*min = 5 N, the upper limit of the cable tension is *t*max = 200 N, *σ*∗ = 1 × 10−6 is the error.

*Xa* = [*θ*1, *θ*2] is the controllable degrees of freedom, and *Xb* = [*θ*3] is the uncontrollable degrees of freedom.

**Figure 3.** 3−6 Static equilibrium workspace of cable-driven parallel robot. (**a**) 3D Static Equilibrium Workspace; (**b**) *xoy* planar view of Static Equilibrium Workspace; (**c**) *xoz* planar view of Static Equilibrium Workspace.

Set the search interval is **S***θ*1 ∈ (65◦, 95◦)&**S***θ*2 ∈ (−40◦, 0◦).

The search step *δθ*1 , *δθ*2 is 1◦, and the pose set to be searched is **Q**.

The static equilibrium working space of the 0−0−2 cable-driven robot is solved by MATLAB programming.

According to Figure 4, the static equilibrium workspace of the 0−0−2 under-constrained cable-driven planar hybrid robot is a curved surface. From Figure 4b, *θ*1*θ*2 planar view of the Static equilibrium workspace, the Static equilibrium points in the pose set are mainly

distributed in the lower right corner area. The uncontrollable degrees of freedom of the static equilibrium point belongs to *θ*3 ∈ (−80◦, −20◦) from Figure 4c,d.

#### **6. Consistent Solution Strategies for Point-to-Point Trajectory Planning**

*6.1. Analysis of the Dynamics of Under-Constrained Systems*

The end poses - *OxO OyO OzO αβγ*T of the parallel robot and Joint coordinates - *θ*1 *θ*2 ··· *θn* T of the planar hybrid robot are written in generalized coordinates *X* = - *x*1 *x*2 ··· *xn* T .

Write Equations (7) and (16) uniformly as:

$$M(X)\ddot{X} - s(X,\dot{X}) - J(X)T = 0\tag{22}$$

where .. *X* is the acceleration in generalized coordinates, *M*(*X*) is the mass matrix, and *J*(*X*)*T* is the vector of the cable tensions.

For the under-constrained system, the controllable degrees of freedom, *Xa*, are *m* dimensional vectors. The uncontrollable degrees of freedom, *Xu*, are *n* − *m* dimensional vectors. As a result, Equation (22) can be written in the form of

$$
\begin{bmatrix} M\_{aa} & M\_{au} \\ M\_{uu} & M\_{uu} \end{bmatrix} \begin{bmatrix} \ddot{X}\_a \\ \ddot{X}\_u \end{bmatrix} - \begin{bmatrix} s\_a \\ s\_u \end{bmatrix} - \begin{bmatrix} I\_a \\ J\_u \end{bmatrix} T = 0 \tag{23}
$$

For Equation (23), when *X*, . *X*, .. *Xa* is known, the acceleration .. *Xu* and the cable force *T* for the system at this moment can be obtained.

By presenting the acceleration terms for the (UDFS) uncontrollable degrees of freedom, we obtain the expressions for .. *Xu*

$$
\ddot{X}\_u = M\_{uu}^{-1} \left( s\_u + f\_u T - M\_{ua} \ddot{X}\_a \right) \tag{24}
$$

Substituting the expression for .. *Xu* into Equation (23), we obtain the expression for the cable extension *T*:

$$T = \left(I\_a - M\_{\rm uu} M\_{\rm uu}^{-1} I\_u\right)^{-1} \left[ \left(M\_{\rm aa} - M\_{\rm au} M\_{\rm uu}^{-1} M\_{\rm uu}\right) \ddot{X}\_a + M\_{\rm uu} M\_{\rm uu}^{-1} s\_u - s\_a\right] \tag{25}$$

#### *6.2. Point-to-Point Trajectory Planning Introducing Parameters*

In the process of point-to-point trajectory planning for under-constrained systems, the start and end points of the motion should be chosen as static equilibrium points and the speed at the start and end points should be zero.

For an under-constrained system with *n* degrees-of-freedom and *m* cables, there are *m* controllable degrees of freedom and *n* − *m* uncontrollable degrees of freedom. The dynamics Equations (22) and (23) shows that during trajectory planning, only controllable degrees of freedom, *Xa*, can be trajectory planned, while uncontrollable degrees of freedom *Xu* are influenced by dynamics [20]. When the controllable degree of freedom reaches the end, an uncontrollable degree of freedom cannot be guaranteed to reach the end with zero velocity. Therefore, it is necessary to choose a suitable trajectory plan for the controllable degrees of freedom to ensure that the uncontrollable degrees of freedom reach the end point with zero velocity.

Assume that the trajectory planned for *Xa* is *xa*(*t*), *t* ∈ [0, *T*], the trajectory of *Xu* to be found is *xu*(*t*), and the set of differential equations is established by Equation (24) as follows.

$$\mathcal{Y} = \begin{bmatrix} \mathbf{x}\_{\mathcal{U}}(t) \\ \dot{\mathbf{x}}\_{\mathcal{U}}(t) \end{bmatrix}$$

$$\dot{\mathbf{y}} = \begin{bmatrix} \dot{\mathbf{x}}\_{\mathcal{U}}(t) \\ \mathcal{M}\_{\text{u}\mathcal{U}}^{-1} (\mathbf{s}\_{\mathcal{U}} + \int\_{\mathcal{U}} T - \mathcal{M}\_{\text{u}\mathcal{U}} \ddot{\mathbf{x}}\_{\mathcal{U}} \end{bmatrix} = f(\mathcal{y}, \mathbf{x}\_{\mathcal{U}}, \dot{\mathbf{x}}\_{\mathcal{U}}, \ddot{\mathbf{x}}\_{\mathcal{U}}) \tag{26}$$

$$\mathbf{y}(0) = \begin{bmatrix} \mathbf{x}\_{\mathcal{U}}(0) \\ 0 \end{bmatrix} := \mathbf{y}\_{0\prime} \ \mathbf{y}(T) = \begin{bmatrix} \mathbf{x}\_{\mathcal{U}}(T) \\ 0 \end{bmatrix} := \mathcal{Y}\tau$$

where *y*0, *yT* are the position and velocity of the non-controllable degrees of freedom at the start and end points, respectively, which are both static equilibrium points with zero velocity.

The optimal trajectory for controllable degrees of freedom is achieved when the trajectory *xa*(*t*) is such that Equation (26) holds.

Since Equation (26) has 2 × (*n* − *m*) variables and its boundary conditions have 4 × (*n* − *m*) constraint. In order to make Equation (26) solvable, 2 × (*n* − *m*) parameter *κ*1, ··· , *κ*2*λ* needs to be added to the planned trajectory *xa*(*t*).

As an example of a conventional polynomial, a trajectory is planned for controllable degrees of freedom *Xa* and parameters are introduced for the planned trajectory.

A conventional polynomial point-to-point trajectory is planned as a straight-line path [21] with a trajectory equation of the form:

$$\mathbf{x}\_a(t) = \mathbf{x}\_a(0) + (\mathbf{x}\_a(T) - \mathbf{x}\_a(0))\mathbf{u}(t) \tag{27}$$

where *xa*(0), *xa*(*T*) is the position of the controllable degrees of freedom at the start and end points, respectively, 0 ≤ *u* ≤ 1.

For a smoothly derivable point-to-point trajectory *xa*(*t*) of order *r*, the law of motion *u*(*t*) is designed as follows:

$$u(t) = \sum\_{i=r+1}^{2r+1} a\_i \left(\frac{t}{T}\right)^i, \ t \in [0, T] \tag{28}$$

where

$$a\_i = \frac{(-1)^{i-r-1}(2r+1)!}{i \cdot r! (i-r-1)! (2r+1-i)!} \tag{29}$$

Introducing parameters to the trajectory of a traditional polynomial programming, the equation of the trajectory is

$$\mathbf{x}\_a(\kappa, t) = \mathbf{x}\_a(0) + (\mathbf{x}\_a(T) - \mathbf{x}\_a(0))\boldsymbol{\mu}(\kappa, t) \tag{30}$$

The law-of-motion function *u*(*κ*, *t*) can be set to

$$u(\kappa, t) = u(\gamma(\kappa, t)) = \sum\_{i=r+1}^{2r+1} a\_i \gamma^i(\kappa, t) \tag{31}$$

where *ai* = (−1) *i*−*r*−1 (2*r*+1)! *i*·*r*!(*i*−*r*−1)!(2*r*+1−*i*)! . In order to maintain *u*(0) = 0, *u*(*T*) = 1, *γ*(κ, *t*). is expressed as follows.

$$\gamma(\kappa, t) = at + \sum\_{i=2}^{2\lambda+1} \kappa\_{i-1} t^i \tag{32}$$

where *α* = 1− 2*λ*+1 ∑ *i*=2 *κi*−1*Ti T* .

Equation (26) is a multivariate marginal differential equation containing parameters that are converted to a multivariate initial differential equation for solution, and the expression of the converted multivariate initial differential equation is as follows.

$$\begin{cases} \dot{y} = f(y, \mathbf{x}\_a(\kappa, t), \dot{\mathbf{x}}\_a(\kappa, t), \ddot{\mathbf{x}}\_a(\kappa, t)) \\ \qquad y(0) = y\_0 \end{cases} \tag{33}$$

Since the solution of the multivariate initial differential equation removes the constraint *y*(*T*) = *yT* at the end point, the result can be solved again by setting up the Newton iteration equation *F*(*κ*) = *y*(*κ*, *T*) − *yT* to ensure that the end point of the solved trajectory is the same as the planned trajectory.

The trajectory parameters *κi* for controlled degree of freedom planning is solved as follows.

#### **7. Simulation Analysis of Trajectory Planning**

#### *7.1. 3*−*6 Under-Constrained Parallel Robot*

Point-to-point trajectory planning was performed for the 3−6 (6 degrees of freedom parallel robots with 3 cables) under-constrained parallel robot with the articulation shown in Table 1, and Table 4 shows the end-effector and trajectory planning parameters.

**Table 4.** End-effector and trajectory planning parameters.

In Table 4, *m* is the mass of the end-effector and *Io* is the inertia tensor of the endeffector in a local coordinate system where the origin of the local coordinate system coincides with the center of mass. *X*0, *XT* is the initial and end pose of the trajectory planning and both are static equilibrium points. Figure 5 shows the initial state and end state for unconstrained parallel robot trajectory planning.

**Figure 5.** Initial state (solid line) and end state (dashed line) for unconstrained parallel robot trajectory planning.

The controllable degrees of freedom are *Xa* = - *OxO* , *O yO* , *O zO* and the uncontrollable degrees of freedom are *Xb* = [*α*, *β*, *γ*]. These are the angles of rotation around the *x*-axis, *y*-axis and *z*-axis respectively.

For the traditional polynomial trajectory planning method, *r* = 3 is taken to ensure that the trajectory is smoothly derivable to the third order, and *ai* is calculated via Equation (29) with the following results.

$$a\_i = \begin{bmatrix} a\_4 \\ a\_5 \\ a\_6 \\ a\_7 \end{bmatrix} = \begin{bmatrix} 35 \\ -84 \\ 70 \\ -20 \end{bmatrix}$$

Substitute *ai* into Equation (27) to find the conventional polynomial trajectory *xa*(*t*), and solve Equation (33) for the multivariate initial value differential equation.

For the polynomial trajectory planning method with the introduction of parameters, take *r* = 3, whose coefficients *ai* are the same as those of a conventional polynomial trajectory. Let the initial value of the parameter vector *κ* be the 6−dimensional zero vector *ζ* = 1 × 10−8, and solve for the parameters using Newton's iterative method to obtain the following results.

Substituting the parameter *κ* into Equation (30), the traditional polynomial trajectory *xa*(*κ*, *t*) is obtained and solved for the multivariate initial differential Equation (33).

The solution results of the traditional polynomial trajectory planning method are compared with those of the polynomial trajectory planning method with parameters, and the comparison results are shown below.

Figure 6 shows the pose of the end-effector in the trajectory. In this case, the direction drops from 0 m to −0.15 m, the direction increases from 0.59 to 0.8 and the direction increases from 1 m to 1.17 m. The orientation angle also changes under the influence of the system dynamics, where *α* rises from −0.43 rad to −0.05 rad and *β* from 0 rad to 0.26 rad (rotated about *z*-axis) and *γ* remains largely stable during the process. The trajectory start and end points of both planning algorithms are consistent with *X*0, *XT* in Table 4 in terms of the pose curves.

Figure 7 shows the velocity curves of the trajectories solved by the two planning algorithms. It can be seen that the difference in velocity between the trajectories planned by the two algorithms is more obvious. Both planning algorithms achieve the boundary condition of zeroing the velocity at the start and end point in the direction of *x*, *y*, *z*. However, for the angular velocity profile of *α*, *β*, *γ*, the traditional trajectory planning algorithm cannot guarantee stationary at the end point, and the trajectory does not zero at *α*, *β*, *γ* of 2.5 s, which cannot guarantee the trajectory boundary condition.

**Figure 6.** *Cont*.

**Figure 6.** Trajectory poses solved under two planning algorithms. (**a**) Comparison of *x*-direction displacement; (**b**) Comparison of *y*-direction displacement; (**c**) Comparison of *z*-direction displacement; (**d**) Comparison of *α* angular; (**e**) Comparison of *β* angular; (**f**) Comparison of *γ* angular.

**Figure 7.** Velocity solved under two planning algorithms. (**a**) Comparison of *x*−direction velocity; (**b**) Comparison of *y*−direction velocity; (**c**) Comparison of *z*−direction velocity; (**d**) Comparison of *α*−angle velocity; (**e**) Comparison of *β*−angle velocity; (**f**) Comparison of *γ*−angle velocity.

#### *7.2. 0*−*0*−*2 Under-Constrained Planar Hybrid Robot*

Point-to-point trajectory planning for an under-constrained cable-driven robot with 0-2-2 (Figure 8). The center of mass of each link is at its geometric center. The coordinates of the articulation points and the link parameters for the cable-driven robot are shown in Tables 2 and 3, while the following Table 5 shows the link inertia and trajectory planning related parameters.

**Figure 8.** Initial state (solid line) and end state (dashed) for under-constrained cable-driven robot in trajectory planning.

In Table 5, *Izi* is the mass product of inertia of the *i* link around the *z*-axis. *X*0, *XT* are the starting and ending points of the trajectory plan, and both are static equilibrium points.

The controllable degrees of freedom are *Xa* = [*θ*1, *θ*2], and the uncontrollable degrees of freedom are *Xb* = [*θ*3].

For the traditional polynomial trajectory planning method, *r* = 3 is taken to ensure that the trajectory is smoothly derivable to the third order, and *ai* is calculated via Equation (29) with the following results.

$$a\_i = \begin{bmatrix} a\_4 \\ a\_5 \\ a\_6 \\ a\_7 \end{bmatrix} = \begin{bmatrix} 35 \\ -84 \\ 70 \\ -20 \end{bmatrix}$$

Substitute *ai* into Equation (27) to obtain the traditional polynomial trajectory *xa*(*t*), and solve Equation (33) for the multivariate initial differential equation.

For the polynomial trajectory planning method with the introduction of parameters, take *r* = 3, whose coefficients *ai* are the same as those of a conventional polynomial trajectory. Let the initial value of the parameter vector *κ* be a 2−dimensional zero vector, *ζ* = 1 × 10−8, and solve for the parameters using Newton's iterative method to obtain the following results.

$$\kappa = \begin{bmatrix} \kappa\_1\\ \kappa\_2 \end{bmatrix} = \begin{bmatrix} 0.58865332\\ -0.29981383 \end{bmatrix}$$

Substituting the parameter *κ* into Equation (30), the conventional polynomial trajectory *xa*(*κ*, *t*) is obtained and solved for the multivariate initial differential Equation (33).

A comparison of the solution results of the traditional polynomial trajectory planning method with those of the polynomial trajectory planning method with the introduction of parameters is shown below.

The joint angle curves for the two planning algorithms are shown in Figure 9. The joint angle *θ*1, *θ*2 starts and ends at the same point in both planning algorithms, while angle *θ*1 decreases from 90◦ to 80◦ and angle *θ*2 changes from −25◦ to −45◦. As can be seen from Figure 7c, in the polynomial trajectory planning method with the introduction of parameters, the start and end points of the *θ*3 curve are the same as *X*0, *XT* in Table 5, whereas the traditional polynomial trajectory planning method does not guarantee this.

**Figure 9.** Trajectory joint angles solved under the two planning algorithms. (**a**) Comparison of *θ*1−angle; (**b**) Comparison of *θ*2−angle; (**c**) Comparison of *θ*3−angle.

The angular velocity of the two planning algorithms is shown in Figure 10. Both planning algorithms reach the boundary condition of zero at the start and end points for the angular velocity of *θ*1, *θ*2. As can be seen in Figure 10c, for the angular velocity profile of *θ*3, stationarity is not guaranteed at the end point using the conventional trajectory planning algorithm, whereas the polynomial trajectory planning method with the introduction of parameters maintains stationarity at the end point.

In summary, for under-constrained parallel or planar hybrid systems, the traditional polynomial trajectory planning algorithm cannot guarantee zero velocity at the end of the trajectory. In contrast, the polynomial trajectory planning algorithm with the introduction of parameters at the start and end points keeps the under-constrained system stable without oscillation and satisfies the trajectory planning requirements.

**Figure 10.** Angular velocity of the solution under the two planning algorithms. (**a**) Comparison of *θ*1−angular velocity; (**b**) Comparison of *θ*2−angle velocity; (**c**) Comparison of *θ*3−angle velocity.

#### **8. Conclusions**

In this paper, the force closure equations and geometric closure equations for incompletely constrained cable traction systems (parallel and planar hybrid) are developed, the equations are solved jointly using a convex optimization solution method with boundary conditions, a static equilibrium inverse kinematic model of the incompletely constrained cable traction system is developed, and the static equilibrium workspace is solved.

Point-to-point trajectory planning algorithms for incompletely constrained cable traction systems are investigated. Traditional purely algebraic methods such as polynomial functions or trigonometric functions for planning the trajectory of the end-effector do not take into account the dynamical model of the system, which can cause oscillations of the system at the start and end points due to insufficient controllable degrees of freedom on the UCR. In this paper, a dynamical model of the UCR is developed, which consider the cable force, external force, controllable position and orientation and uncontrollable position and orientation of the end-effector, etc. Based on this model, the traditional polynomial-based point-to-point trajectory planning algorithm is improved by adding 2 × (*n* − *m*) parameters to the kinematic law function *u*(*t*). The constraints of the dynamics model are incorporated into the trajectory planning process to achieve point-to-point trajectory planning for the UCR. The results show that the trajectory of the improved algorithm is smooth and derivable, and the end-effector is stationary and stable without oscillation at the start and end points, which proves the effectiveness of the optimized algorithm.

**Author Contributions:** Conceptualization, Q.D., Q.Z. and T.W.; methodology, Q.D.; software, Q.Z.; validation, Q.D. and Q.Z.; writing—original draft preparation, Q.Z.; writing—review and editing, Q.D. and Q.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Key Research and Development Program of Shaanxi (Program No. 2022GY-314).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**

## *Article* **Singularity Analysis and Geometric Optimization of a 6-DOF Parallel Robot for SILS**

**Doina Pisla 1, Iosif Birlescu 1,\*, Nicolae Crisan 2, Alexandru Pusca 1, Iulia Andras 2, Paul Tucan 1, Corina Radu 3, Bogdan Gherman 1 and Calin Vaida 1**

	- 400000 Cluj-Napoca, Romania

**Abstract:** The paper presents the singularity analysis and the geometric optimization of a 6-DOF (Degrees of Freedom) parallel robot for SILS (Single-Incision Laparoscopic Surgery). Based on a defined set of input/output constraint equations, the singularities of the parallel robotic system are determined and geometrically interpreted. Then, the geometric parameters (e.g., the lengths of the mechanism links) for the 6-DOF parallel robot for SILS are optimized such that the robotic system complies with an operational workspace defined in correlation with the SILS task. A numerical analysis of the singularities showed that the operational workspace is singularity free. Furthermore, numerical simulations validate the parallel robot for the next developing stages (e.g., designing and prototyping stages).

**Keywords:** parallel robot; singularity analysis; geometric optimization; single-incision laparoscopic surgery

#### **1. Introduction**

Single-incision laparoscopic surgery (SILS) is a type of minimally invasive surgery (MIS) where the surgeon uses a special multi-lumen trocar (a trocar with multiple insertion points) inserted through a single incision in the patient body for all instruments required for surgery. SILS was described in the gynecological literature as early as 1969, when Wheesess reported the first 4.000 cases of tubal ligation [1]. As the literature suggests, the main advantages of SILS are: (i) the short recovery time for the patient, (ii) lesser degree of pain (when compared to classical MIS), (iii) and better cosmetics (since it uses a single access port) [2]. However, the challenges of SILS (performed manually) are caused by its poor ergonomics [3]. As history shows, robotic systems for surgery were developed to circumvent various limitations of the classical, human-performed interventions. MIS is one example in which robots had a significant impact [4,5]. Multiple approaches have been proposed and developed for robotic systems for MIS, namely: multi-arm robotic systems [4], flexible robotic systems [4], laparoscope holders [6], etc. The analysis of the intraoperative data from robotic surgery has enabled the addition of new technologies, one of them referring to automated gesture recognition, which assists the surgeon during delicate procedures such as suturing. For this task, multiple intelligent algorithms were developed using supervised, unsupervised, or semi-supervised learning paradigms [7]. This is an important step toward autonomous robotic-assisted surgery, where the robot uses a confidence-based shared control strategy to perform certain tasks under supervision but without the involvement of the human operator [8]. The advancements in computer science, which are now providing more reliable intelligent algorithms and the increased availability of large quantities of information (big data)—pre-operative digital information,

**Citation:** Pisla, D.; Birlescu, I.; Crisan, N.; Pusca, A.; Andras, I.; Tucan, P.; Radu, C.; Gherman, B.; Vaida, C. Singularity Analysis and Geometric Optimization of a 6-DOF Parallel Robot for SILS. *Machines* **2022**, *10*, 764. https://doi.org/10.3390/ machines10090764

Academic Editor: Zhufeng Shao

Received: 29 July 2022 Accepted: 29 August 2022 Published: 2 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

instruments motion for repetitive tasks, recorded surgeries, patient evolution in the short and long term—are indicators and supporters of the next generation of devices, intelligent surgical robots [9].

The year 2018 marks the starting point for surgical robots dedicated to SILS due to the first use of the robotic da Vinci SP system in SILS surgery [4,10,11]. The da Vinci SP robot received FDA approval in 2018, becoming a benchmark in SILS robotic surgery, with a serial architecture capable of manipulating and orienting active instruments and the laparoscopic camera by using a single port through which these instruments are inserted. The system has an articulated laparoscopic camera and instruments, thus offering much better visibility and dexterity in the operating field compared to classical instruments [12]. The drawbacks are the limited forces (for tissue manipulation) due to the multiple bends of the instruments. Another robot dedicated to SILS surgery that has received FDA approval is the Senhance robot, which in contrast to da Vinci SP is a multi-arm system with three independent serial manipulators capable of orienting and manipulating surgical instruments [13]. Furthermore, there are various robotic systems for SILS which have not yet received FDA approval, such as SORT, SPIDER, MASTER, Virtual Incision, and so on [4].

Since a robotic system for SILS must comply with complex design specifications (e.g., footprint in the operating room, occupied volume, quick removal in case of an emergency), various optimizations are achieved in the early design stages. In a general interpretation, optimization aims to find the best solution for a given problem based on the imposed restrictions. In robotics, the optimization process often faces conflicting criteria such as speed versus energy efficiency, accuracy versus stability, or geometric dimensions versus working volume. When dealing with multi-objective optimization problems, the results are often given as a set of Pareto optimal solutions which must be further interpreted. In [14], the authors present an optimal dimensional synthesis of a parallel mechanism using two objective functions: finding the smallest dimensions for the geometric links of the robot for a given workspace, and the second aims to ensure the best overall dexterity within this workspace. As the two functions have conflicting criteria, it was shown that favoring the first function can lead to a robot with low dexterity while the second leads to a bulky robot. A numerical optimization methodology is presented in [15] to achieve an optimal design synthesis for a planar parallel manipulator for a prescribed dexterous workspace, using the condition number (an index that describes the dexterity of a robot). Merlet emphasizes in [16] the importance of the proper use of the condition number and global conditioning index (GCI) in the optimal design of a robot, illustrating their limitations (namely, when they are applied to robots that have both translation and orientation motions). GCI was first introduced by Gosselin and Angeles in [17], and it was one of the first approaches to describe the global performance of a manipulator, also demonstrating that in some cases, GCI can provide conflicting results in workspace optimization problems. In [18], the authors successfully applied the GCI and condition number to a 2-DOF parallel mechanism enabling them to obtain a generalized characterization of the manipulator with respect to its operational workspace. In [19], the authors used the condition number to optimize Orthoglide, a 3-DOF parallel robot with only translational motions used in milling applications. The dimensional synthesis of the 3- DOF Delta parallel robot for a prescribed workspace is presented in [20], where the authors define several optimization objectives solved using the Lagrange multipliers method, demonstrating that existing industrial robots could have almost half their size for the same given workspace. Furthermore, for a Delta-like parallel robot, in [21], the authors achieved an optimum design using multi-objective optimization algorithms, i.e., Pareto-optimization. A methodology to achieve an optimal design for a 6-DOF parallel manipulator having as objective its accuracy is presented in [22], where the authors are using the distribution of the condition number to determine the best solutions. The Stewart–Gough platform was analyzed in detail in [23], where the authors used a multi-objective optimization algorithm (NSGA-II) to obtain the optimal geometric parameters and leg stroke lengths, demonstrating that this approach is more efficient than the classical numerical methods. The workspace maximization of the Delta robot using a geometrical technique implemented in CATIA is shown in [24], where the use of the Simulated Annealing Algorithm provided promising results. The process of optimization is also found in the scientific literature concerning the development of medical robots. In [25], the authors optimized a 3-DOF serial robot considering the robot workspace and the surgical instrument insertion points. Other authors [26] used optimization processes for the static balancing of a surgical robot.

A novel approach for robotic-assisted SILS was proposed in [27], which is based on hybrid robotic systems with three major components: (1) a 6-DOF parallel robot that simultaneously guides (2) two interconnected orientation platforms (mounted on the mobile platform of the 6-DOF parallel robot) together with a laparoscopic camera, and (3) two active surgical instruments (mounted in the orientation platforms) to perform the SILS task. The kinematics of the surgical instruments was studied in [28], and in [29], the input–output equations for the orientation platform were provided. Furthermore, a recent work [27] achieved the kinematic modelling of the 6-DOF parallel robot.

The present paper illustrates the results regarding the further development of the SILS robotic system, namely: (i) a singularity analysis for the 6-DOF parallel robot presented in [27] (correlated with the operational workspace of the robot) and (ii) a geometric optimization of the 6-DOF parallel robot, to reduce its operating room footprint while maintaining the capabilities of performing the SILS task in a defined operational workspace. The singularity analysis was achieved using the vanishing conditions of the Jacobian matrices A and B. The input singularities are straightforward to interpret for the 6-DOF parallel robot, and for the output singularities, a geometric interpretation was achieved based on the characteristic tetrahedron [30]. For the geometric optimization, the Multi-Objective Genetic Algorithm implemented in Matlab gamultiobj function [31] is used to minimize the dimensions of the robot for a defined operational workspace. The GCI index is used in the optimization process only as guidance and not as a decision parameter due to its drawbacks when used on robots that have both translation and orientation motions [16,17].

The paper is structured as follows: Section 2 defines the parallel robotic system with all its modules and presents the inverse geometric model for the 6-DOF parallel robot. Section 3 defines the robotic system's operational workspace to create a clear description of the robotic system task. Section 4 presents a singularity analysis for the general 6-DOF parallel robot. Section 5 presents the proposed optimization algorithm for the 6-DOF parallel robot and the general optimization criteria. Section 6 presents the optimized version of the 6-DOF parallel robot, showing that the operational workspace is singularity free. Furthermore, numerical simulations are provided to validate the optimized robotic system. The conclusions are presented in Section 7.

#### **2. An Innovative Parallel Robotic System for SILS**

SILS requires the simultaneous independent manipulation of two surgical instruments while the operating field is viewed through a laparoscopic camera. The position of the laparoscope is not fixed, and it can change depending on the target tissue location. Furthermore, both the surgical instruments and the laparoscope must be manipulated with respect to Remote Center of Motion (RCM) points, i.e., the insertion points of the instrument within the trocar. Lastly, the position of the insertion points (or the trocar) on the patient's abdomen may vary based on various medical factors (e.g., the position of the resected tissue, areas that must be avoided, etc.) [3].

To comply with the required design parameters, an innovative parallel robotic system for SILS was proposed with the following components:

3. **Two surgical instruments** (described in [28])—each mounted in one orientation platform. Each surgical instrument has a serial architecture with 4-DOF: 1-DOF for the rotation about its longitudinal axis (the third rotation for the surgical instrument), 1- DOF for the (articulated) bending, and 2-DOF for the gripper (each jaw of the gripper is actuated separately enabling grabbing and gripper turning) (Figure 1c).

**Figure 1.** The initial modules for the SILS robotic system: (**a**) the 6-DOF parallel robot; (**b**) mobile platform containing the laparoscope and two orientation platforms; (**c**) surgical instrument.

With respect to a defined robotic-assisted medical protocol [3], the proposed robotic system must comply with the following stepwise procedure:

The kinematic scheme of the 6-DOF parallel robot is presented in Figure 2. Although the mechanism topology is the same as in [27], in this work, the kinematic chains are not considered identical (removing this design constraint may improve the final solution for the robot architecture). In [27], the links *l*2, *l*3, and *l*4 had the same values.

**Figure 2.** Kinematic scheme of the 6-DOF parallel robot for SILS.

The 6-DOF parallel robot mechanism topology is as follows [27]:

#### *Inverse Geometric Modelling for the 6-DOF Parallel Robot*

The inverse geometric model will provide constraint equations which, in turn, will be used as objective functions for the optimization algorithm.

For the inverse geometric model, the inputs are the Cartesian coordinates of point *E* [*XE, YE, ZE*] and the orientation angles *ψ, θ, ϕ*, whereas the outputs are the active joint generalized coordinates *qi* (*i* = 1 ... 6). With respect to the *ZYX* Euler convention, the coordinates of the passive spherical joints are: