Multi-Representations and Multi-Constraints Approach to the Numerical Synthesis of Serial Kinematic Structures of Manipulators

This paper presents a set of algorithms for the synthesis of kinematic structures of serial manipulators using multiple constraint formulation and provides a performance comparison of different kinematic representations, the Denavit-Hartenberg notation, the Product of Exponentials (screws), and Roll-Pitch-Yaw angles with translation parameters. Synthesis is performed for five given tasks, and both revolute and prismatic joints can be synthesized. Two different non-linear programming optimization algorithms were used to support the findings. The results are compared and discussed. Data show that the choice of the constraint design method has a significant impact on the success rate of optimization convergence. The choice of representation has a lower impact on convergence, but there are differences in the optimization time and the length of the designed manipulators. Furthermore, the best results are obtained when multiple methodologies are used in combination. An arbitrary manipulator was designed and assembled based on a trajectory in the collision environment to demonstrate the advantages of the proposed methodology. The input/output data and synthesis methodology algorithms are provided through an open repository.

[15] combined their analytical approach with numerical op-80 timization. Today, researchers tend to use optimization algo-81 rithms to synthesize manipulators for complex tasks. Mathe-82 matical optimization is a set of numerical methods that search 83 for an optimal solution that meets a given objective function 84 while satisfying constraints, if specified. Mostly, global or 85 local minimum search algorithms are applied. Among global 86 optimization algorithms, in the case of synthesis of a robot 87 manipulator, Genetic Algorithms have been applied in vari-88 ous studies. The comparison of straight, rounded, and curved 89 links was studied by Pastor et al. [16], using Schunk actuator 90 modules as joints. The design parameters vary depending on 91 the type of link, without directly specifying any kinematic 92 representation. Valsamos and Katrantzis et al. [17], [18] 93 also used Schunk actuators with pseudo-joints between them 94 and generated the optimal kinematic structure based on the 95 manipulability of the manipulator. Another global optimum 96 search algorithm is the Simulated Annealing implemented in 97 [19], where they composed a manipulator from predefined 98 links and joints models. Singh et al. [20] applied Binary 99 Search Algorithm to synthesize DH parameters to assembly 100 a modular manipulator. 101 The studies mentioned above were based on the synthesis 102 of predefined motor models, and their approach is limited 103 to those actuators; that is, the output kinematic structures 104 are not purely arbitrary. The synthesis of arbitrary structures 105 was investigated by Patel and Sobh [21] who searched with 106 Simulated Annealing Algorithm the optimal set of  Hartenberg (DH) parameters for a given goal and then tested 108 the inverse kinematics with Particle Swarm Optimization. 109 Sigla et al. [22] obtained DH parameters of an arbitrary 110 manipulator using the Augmented Lagrangian method. 111 The main disadvantage of global optimum search algo-112 rithms is that, for complex tasks, they can search for a 113 solution within minutes, hours, or days, even on today's pow-114 erful computers. Therefore, local optimization algorithms 115 can serve better, especially at the beginning of the custom 116 manipulator design process. Dogra et al. [23] utilize a non-117 linear programming method based on the minimization of 118 joint torques, which means that it is a dynamical synthesis 119 and the objective is not to reach given poses, but to fulfill cri-120 teria based on energy consumption. The design parameters do 121 not follow any standard kinematic representation, since they 122 include dynamics parameters. Whitman et al. [24] search 123 for the optimal design of a robotic arm with an objective 124 function that minimizes the length of the path in the joint 125 space. The design parameters are presented vaguely without 126 any specified standardization as a matter of choice and can in-127 clude "link lengths, twists about link axes, base locations, or 128 other offsets". Shirafuji and Ota [25] implemented a synthesis 129 method based on differential inverse kinematics, where they 130 avoided dual optimization of inverse kinematics and robot 131 design parameters. This contrasts with the majority of applied 132 methods, including the approach presented in this paper. 133 However, one of the disadvantages is that the orientation of 134 the goal pose is not taken into account. The optimization 135 constraints are based on Jacobian, and the outputs are screw 136 coordinates. 137 Based on our literature review, two major questions are 138 raised and will be addressed in this study. At first, there 139 are multiple robot kinematic representations; however, the 140 2 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10  to that, we added a method that does not apply the norm 167 of the transformation matrix but its elements, a method that 168 uses the norm of the dual quaternions, and we also applied a 169 new method given as a metric in a 12-dimensional Euclidean  The methodology is presented in Section III. The results 189 are summarized and discussed in Sections IV and V. We process; that is, the same input provides the same output.

194
In Section VI we present the proof-of-concept of the  revolute (3R) mechanism that was built based on a specified 196 collision environment. In this section, one can also find a per-197 partes optimization process that leads to the final manipulator 198 design. This contrasts with other studies that synthesize 199 kinematics, dynamics, and collision avoidance at the same 200 time. A simple method for avoiding joint collisions with the 201 environment is presented.

203
This section is divided into three parts. In the beginning, the 204 mathematical background is briefly introduced. The descrip-205 tion of the differences between the optimization of DH, PoE, 206 and RPYXYZ follows. The last part describes the methods 207 for creating constraints.

208
For optimization itself, the MATLAB TM function fmincon 209 [27] was chosen, which is a non-linear programming solver. 210 It searches for the local minimum of an objective function 211 within defined constraints. Beside others, it supports two 212 optimization algorithms that were suitable for the synthesis 213 problem -interior-point (IP) [28] and sequential quadratic 214 programming (SQP) [29].

215
In this study, 3780 optimizations of the kinematic struc-216 tures were performed, as shown in Table 1. The number 217 of representations is three, and the number of constraint 218 types is the combination of six methods, i.e., sixty-three 219 possibilities from using every method individually to using 220 all six combined. Five different paths were chosen; they are 221 shown in Figure 1. For every path, the placement of a base 222 of the robot was restricted in two ways -in the world frame 223 (that is, the position of the base was not optimized) and in 224 the interval [−0.3, 0.3] [m] in every direction around the 225 world frame. Those two options were chosen to demonstrate 226 how the possibility of adjusting the base of a robot during 227 optimization can influence its convergence. One can also see 228 it as a different task (a new given path) when the searched 229 kinematic structure can be completely different. For each optimization run, the synthesis algorithm requires 231 specifying these inputs by a user:  The open-source Matlab toolboxes prepared by Lynch [9] 244 and Corke [30] were used for the calculations presented.

245
The MathWorks Robotics System Toolbox was used for 246 visualization, inspection, and verification.
Furthermore, the one-to-one mapping of SE(3) in the 256 projective space P(R) 7 was used over the eight-dimensional 257 vector space of dual quaternions [31]. Dual quaternion p consists of which is the quaternion representing orientation, and which is the dual part representing the position in space. 261 Together, they form homogeneous coordinates (a 0 : a 1 : a 2 : 262 a 3 : c 0 : c 1 : c 2 : c 3 ) in P(R) 7 , where the elements of 263 SE(3) can be represented by points. One of the benefits of 264 this representation is that both position and orientation can 265 be expressed in a single vector.

266
Joints can be kinematically expressed using homogeneous 267 transforms, or as normalized twists called the screw axis S, 268 which is 6 × 1 vector The vector ⃗ ω is the direction of the rotational axis and 270 corresponds to the vector ⃗ a of a transformation matrix. A 271 coordinate q is any point on the axis. Their cross product 272 ⃗ v = −⃗ ω × q forms a moment of the axis about origin. Note 273 the difference with the Plücker coordinates where ⃗ v ′ = ⃗ ω ×q. 274 The forward kinematics of a manipulator is calculated using 275 the product of exponentials formula, as shown in 6. The where, in the case of DH convention, is obtained by multiplying the rotation matrix R z (θ i ) The displacement calculated with the RPYXYZ parame-293 ters is where T i is the translation matrix between the origins 295 of the two frames and R xi , R yi , R zi are rotations around In the options of the fmincon solver, the maximum itera-321 tion value was increased to 300. The number of maximum 322 function evaluations was set to 30, 000. The Step Tolerance 323 was set to 10 −300 . The other function options were kept as 324 default values. The fmincon requires an initial estimation of 325 the optimized parameters. For a local optimum search, this 326 first guess may have a major impact on the outcome. To avoid 327 any randomness, we applied the methodology D described in 328 our previous paper [32], where four methods are presented to 329 create kinematic structures represented in the DH parameters 330 using geometrical analysis between a given pose and the base 331 of the manipulator. The methodology D places the joints 332 according to the DH convention on a Bézier curve. The last 333 joint is placed out of the curve but in a way that the forward 334 kinematics of the serial chain can reach the chosen pose, in 335 this study, the most distant pose on the path. The distance 336 is measured from the center of the given interval where the 337 base of the robot can be placed. The fmincon solver also 338 requires the specification of the bounds. These limits were set 339 as presented in Table 2. The length and moment parameters 340 are determined again in relation to the most distant pose of 341 a given path and the norm of its position vector ⃗ t. The dual 342 quaternions are in projective space; therefore, the bound is 343 set to infinity.

344
The base of a robot is represented as a dual quaternion p b 345 with the orientation part set to the identity (p b0 = 1; a 0 = 346 1, a 1 = 0, a 2 = 0, a 3 = 0), that is, no rotation is allowed, and 347 only the dual part p b1 (base translation) is being optimized. 348 The end-effector is represented as dual quaternion p e without 349 any optimization restrictions. These are in total the first 350 twelve optimization parameters that enter the optimization 351 process: 352 n be = n b + n e = 4 + 8 = 12 (10) Note that if no base interval is given, then the robot base 353 frame displacement is not optimized, but set to the input 354 value or the world frame.

355
The following subsections describe how the algorithm 356 design varies according to the chosen robot kinematic rep-357 resentation. 358

1) Denavit-Hartenberg parameters 359
The objective function is to minimize the length parameters 360 given by the DH convention, that is, the absolute distances 361 VOLUME 4, 2016 a i and d i between the axes z i−1 to z i and x i−1 to x i of the 362 joints, as shown in (11).
where n j is the number of joints and ⃗ t e is the 4 × 1 vector The objective function for the screw representation is de-382 signed as follows.
where ⃗ m is 4 × 1 vector of the dual part (position) of the In addition to that, every pair of vectors ⃗ ω i and ⃗ v i preserves 399 the perpendicularity between these vectors and their dot 400 product shall be equal to zero.
Therefore, equations (13) and (14) are also added to the set 402 of optimization constraints. The objective function for the RPYXYZ parameters is set where ⃗ t e is again the dual part of the end-effector displace-406 ment in relation to the last joint frame, ⃗ t i is the vector of the 407 translational coordinates XYZ between consecutive joints.

409
Only the type of equality constraints of the fmincon solver 410 are applied in the presented methodology, and the inequal-411 ity constraint vector remains empty. The algorithm tries to 412 minimize the objective function (provide as short linkages as 413 possible) which is expected to not be zero in the end, but on 414 the other hand, it has to satisfy the given constraints, i.e., the 415 constraint violation is not allowed and the sum of constraint 416 errors has to be minimized to zero.

417
Let us define the sum of errors E -the constraint violation 418 value as the sum of all constraints equality equations ceq.
which is equivalent to the difference between the forward 420 kinematics of the synthesized arm J n,k (θ n,k ) and the n given 421 poses P n,k . The number of constraining equations k depends 422 on the chosen method of constraint creation and, because 423 they can be combined, also on their combination.

424
For each given pose P 1..n , the constraints are calculated 425 in every iteration and based on this output, the kinematic 426 structure design parameters are altered and parsed in the next 427 iteration.
there is a set ceq A of twelve constraint equations defined as where R i represents i = 1..9 th element of rotational 433 part of matrix J or P, ⃗ t represents x, y, z elements of the 434 translational vector of homogeneous matrices J or P. w is 435 the weight that compensates the ratio between meters and 436 radians.
This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.
where the vectors ⃗ n, ⃗ o, ⃗ a represent x, y, z axis vectors of where p J(θ) and p P are dual quaternion representation of which is the norm between those two dual quaternions.

FIGURE 2.
Visualization of an optimized kinematic structure performing Path 5; the base and the end-effector frames are pink; the joint frames have RGB color, joint rotational axis is blue.

468
This section presents an analysis of the output data from 469 the optimization process. The input data (paths), multi-470 representations and multi-constraints synthesis algorithm, 471 and the results are available on public repository [34]. Based 472 on the analysis of obtained data, the performance of the 473 representations is similar; however, the choice of constraint 474 design seems to have a major impact on the overall synthesis 475 problem. An example of an output kinematic structure from 476 a single optimization run out of those 1890 is visualized in 477 Figure 2.

478
The convergence success rate is presented. It shows how 479 often the algorithm provides an output with the constraint 480 violation value E < 10 −5 , i.e. the structure is able to 481 precisely reach all poses on a given path; see again Equation 482 (16) for details. Another parameter that is compared is the 483 length of the optimized manipulator, calculated as the sum 484 of the distances between consecutive joint axes given by 485 the nearest points on those axes. Instead of the number of 486 iterations, which did not differ much for the same inputs, the 487 optimization time is provided. The length of the manipulator 488 and the optimization time are also presented only for the 489 outputs with E < 10 −5 .  Table 3 shows how the complexity of a given path can influ-492 ence the output, i.e. when the given task is more complicated 493 from the synthesis point of view, the convergence success rate 494 drops.

495
In Figure 3(a), the performance of robot kinematic repre-496 sentations can be observed. In general, there are no unex-497 pected differences between them. For every representation, 498 630 optimizations were performed. As already mentioned 499 in the Introduction part, the RPYXYZ was not used in any 500 VOLUME 4, 2016

Manipulator length [m]
(c) Performance of constraint combination methods for all paths; convergence success rate, optimization time, manipulator length. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.

532
The output data in Figure 3( Figure 4, where 546 the convergence success rate is shown only for complicated 547 (from the synthesis point of view) Paths 2, 4, and 5. Although 548 some methods keep their success rate similar regardless of 549 what task is given, many other methods do not achieve even 550 a single successful run.

551
Another interesting observation is that the combination of 552 constraint methods neither extends the optimization time nor 553 seems to have an impact on the length of the manipulator. 554 However, based on the data obtained, the combination of 555 more than three methods does not produce a higher success 556 rate.

558
In Figure 5, the same data are provided for the SQP algo-559 rithm. The average manipulator length is similar, the conver-560 gence success rate is slightly lowered, and the optimization 561 time is extended a lot. Therefore, the IP algorithm provides 562 better results than SQP.

563
The most interesting comparison provides the plot of the 564 convergence success rate for the constraint design methods. 565 The performance of particular methods A to F is reduced. 566 The most successful combinations of AC, AD and ACD also 567 decreased, but not so much. This leads to the suggestion that 568 the combination of AC, AD, or ACD methods should be 569 used as the first choice in future research on the problem of 570 synthesis of kinematic structures.

571
The serial manipulator synthesis problem using numerical 572 optimization is a very challenging and time-consuming task 573 where many variables can have an impact on the outcome. 574 Therefore, we do not want to claim that constraint design 575 method A and its combination will always be the best choice 576 VOLUME 4,2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and

Manipulator length [m]
(c) Performance of constraint combination methods for all paths; convergence success rate, optimization time, manipulator length. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3186098 for the synthesis of robot kinematic structures. It can be a 577 coincidence that these methods suffer less in terms of singu-578 larities in combination with the chosen optimization solver 579 fmincon and the related IP and SQP algorithms. This shall be 580 verified in future studies; however, we want to highlight this 581 problem because the constraints design methodology has not 582 been studied before and seems to be of particular importance. the convergence success rate. Note that it was still achieved 631 by avoiding any randomness in the optimization process. 632 This fast iterative procedure may speed up the overall time of 633 the synthesis and custom manipulator design process. While 634 the DH convention and the RPYXYZ parameters share the 635 idea of partial transformations between joints,the PoE uses a 636 completely different computational approach for direct and 637 inverse kinematics problems, so switching between these 638 representations during synthesis, especially in local minima 639 search algorithms, can help overcome a local minimum to get 640 closer to the global minimum.

641
There are other variables that need to be investigated in 642 more detail. We chose the weight w = 0.1 to compensate 643 the relation between meters and radians based on our expe-644 rience with overall performance. Definitely, every constraint 645 method would perform differently if an exact value was de-646 termined for each separately on the basis of its performance. 647 However, the value of w can also vary based on the number 648 of poses of a given path and the overall size of the trajectory 649 (and output manipulator). This requires another broad study. 650 On the other hand, even though in some presented methods 651 the translantion may seem preferred over the orientation 652 with this weighing, we strived to reduce its impact using a 653 very strict convergence threshold (sum of positioning errors) 654 E < 10 −5 .

655
So far, six constraint design methods have been examined. 656 These can be extended to other ones, e.g. the combination of a 657 quaternion for rotation and XYZ coordinates for translation, 658 their norms, or even other metrics such as Method C. Also 659 there is a possibility to test performance of other norm types 660 than 2-norm, for example the Frobenius norm.

662
This chapter describes how the presented methodology can 663 be applied in the design of serial robotic manipulators. After 664 the optimal kinematic structure is synthesized, there are a 665 few challenges related to dynamics and collision avoidance 666 problems. In contrast to other studies [16], [19], [20], [22]-667 [24] that try to synthesize kinematics along with dynamics 668 and/or collision avoidance at the same time, we propose a 669 per-partes optimization approach. This means dividing the 670 optimization process into steps: 671 1) Early synthesis of kinematic structure to obtain better 672 initial estimation of the optimized values. 673 2) Synthesis of kinematic structure using the output from 674 the previous optimization as initial estimation, chang-675 ing representations and/or constraint design methods to 676 achieve optimal result.  (b) DH parameters and constraint method A as initial estimation.
Convergence success [%] (d) DH parameters and constraint method A as initial estimation.  in other studies, for example, in the case of links [37] and in 697 the case of choosing the right actuators [38].

699
The fmincon function was implemented again with a similar 700 goal as in the previous ones, the objective being to keep the 701 links as short as possible: where || ⃗ t i || is the distance between two consecutive joints 703 in RPYXYZ representation. The start index is i = 3, which 704 is the displacement between the first and second joints in this 705 representation. The displacement between the world frame 706 and the robot base frame has i = 1, the displacement between 707 the base frame and the first joints with i = 2 is expected to 708 be in a position without any collision.

709
The physical position of every joint i represented by a 710 homogeneous matrix J i can be expressed as a point r i on 711 its rotation axis using the parametric equation of a line: with parameter s i . The vector ⃗ t represents the position and 713 ⃗ a the vector of the z axis of J i . Equation (33) therefore forms 714 a set of constraint equations. Another constraint set deals 715 with collisions. The environment is represented by spheres 716 of diameter d s and the joints by spheres of diameter d j .
where k is the number of pairs to be compared -every 718 joint sphere with every environment sphere. || ⃗ t j − ⃗ t s || is 719 the distance between those two spheres in the world frame. 720 This methodology for creating ceq k constraints can be re-721 placed by some more advanced collision check methods [39]. 722 Constraints defined by Equations (33) and (34) where B is the position of the base frame of the robot, 743 M is the displacement between the base and the end-effector 744 frame, and S 1..3 are rotational screws.

745
When reaching pose 2, the structure collides with one 746 of the walls of the environment as shown in Figure 8(b). 747 The methodology described in the previous subsection was 748 applied to find the collision-free position of the joints. The 749 joint position optimization result in a simplified collision 750 environment (spheres only) is visualized in Figure 8(c).

751
The final robot design is visualized in Figure 9. The links 752 between the joints were manually modeled and 3D printed. 753 The motors Dynamixel XH430-V350-R were chosen for 754 assembly in an experiment to perform the given trajectory. 755 A video of the robot performing the given task is attached to 756 this paper and may be found in the supplementary material.  ture. Using this multi-representation and multi-constraint 800 approach, it is possible to apply a relatively good result 801 from one optimization, where, for example, DH parameters 802 were searched for, convert them to PoE, and start over with 803 a kinematic structure that is already close to the optimal 804 result, but calculating with a very different approach -matrix 805 exponentials and screw coordinates instead of partial trans-806 formations. This increases the convergence success rate and 807 decreases the length of the output manipulator along with the 808 time needed for optimization.

809
To verify the results presented, an arbitrary manipulator 810 was built, and an experiment was carried out along with the 811 introduction of the split optimization process. The complex 812 task of finding the optimal manipulator design was divided 813 into five steps, where the task of finding the kinematic struc-814 ture, collision-free performance on a given trajectory, and 815 dynamics analysis was either demonstrated in this paper or 816 referenced in existing studies.

818
The cooperation between the authors was initiated thanks 819 to the Aktion Osterreich-Tschechien, AOCZ Semesterstipen-820 dien grant, financed by Federal Ministry of Education, Sci-821 ence and Research (BMBWF) and organized by Austrian 822 Agency for International Cooperation in Education and Re-823 search (OeAD-GmbH).