GlobDesOpt: A Global Optimization Framework for Optimal Robot Manipulator Design

Robot design is a major component in robotics, as it allows building robots capable of performing properly in given tasks. However, designing a robot with multiple types of parameters and constraints and defining an optimization function analytically for the robot design problem may be intractable or even impossible. Therefore black-box optimization approaches are generally preferred. In this work we propose GlobDesOpt, a simple-to-use open-source optimization framework for robot design based on global optimization methods. The framework allows selecting various design parameters and optimizing for both single and dual-arm robots. The functionalities of the framework are shown here to optimally design a dual-arm surgical robot, comparing the different two optimization strategies.


I. INTRODUCTION
Properly designing robots is a time-consuming and challenging task, which requires great expertise and knowledge on the field of application of the system. A good design should ensure good performance regardless of the task that it is given, otherwise different designs should be needed for different tasks. Thus, robotic design optimization with consideration of different types of parameters and constraints is of great significance and practicability.
Several measures can be employed to define an optimal design, such as the performance in executing a desired control task, accuracy in path tracking, and robot dexterity. Solvers like [1] using standard optimization techniques can be used for robot kinematic calibration and to improve the robot model, however, analytically formulating an optimization function may be challenging or even impossible, especially for highly complex robotic structures. For these reasons, researchers have focused on black-box optimization methods, which do not require a close form of the optimization function, nor knowledge about its derivatives. For instance, Hassan et al. [2] employed Genetic Algorithm (GA) to optimize the design of a gripper; in [3] the optimal design of a 7DOF robot is studied using GA and considering different The associate editor coordinating the review of this manuscript and approving it for publication was Luca Cassano. kinematic and dynamic measures; Particle Swarm Optimization (PSO) was used in [4] to optimally design a cable-driven robot and [5] employed Adaptive Simulated Annealing for the design of a concentric tube robot. Most of the works focus on the design of robots with just one arm. However, in fields like minimally invasive surgery, surgical robots are usually employed in a dual-arm or even multiple-arms configuration [6]- [8], both at the master and slave side. Optimally designing dual-arm robots is more challenging, as the interaction between the two arms should be considered [9]. Lee et al. [10] propose a dual arm manipulability measure specifically designed for multi-arm systems, and Torabi et al. [11] extended it to master-slave teloperation. Only a few works in the literature have focused on the optimal design of dual-arm robots, such as [12] where PSO is used to design a dual-arm concentric-tube robot, or in [13] GA for a bi-manual articulated robot. However, both works focus on task specific applications.
Recently, in the field of machine learning, Bayesian Optimization (BO) has been gaining a lot of popularity and proved very efficient in solving black-box optimization problems. It is considered the state-of-the-art machine learning framework for optimization tasks in terms of data efficiency, and has been successfully applied in engineering, machine learning, and design [14], [15]. BO has also been widely used in robot control to optimize control parameters and policies [16]- [18] as it makes efficient use of data and experiments by learning a probabilistic surrogate model of the function to optimize, which is used for finding optimal parameters. By exploiting the learned model, BO requires fewer interactions (i.e., evaluations of the true objective function) than other optimization methods [16].
Lately, in [19] a method for optimizing the design of a single surgical robotic arm was proposed, yet the approach results to be little flexible, as it is very application-specific and time consuming. To the best of the authors' knowledge, no software is available online currently to simplify the robotic optimizing design process. We have therefore developed GlobDesOpt (available at https:// github.com/cursi36/GlobDesOpt) a versatile and flexible Matlab package for optimal robot design, based on three different solvers (PSO, GA, and BO).
The contribution of this manuscript is two-fold: • providing GlobDesOpt, a simple-to-use open-source Matlab package for optimal robot design; • proposing two optimization strategies which can be applied to both single and dual-arm robots, with the dual-arm optimization including the dual-arm dexterity measure and a safety measure to maximize the dual-arm robot's dexterity while reducing the risk of collisions between the two arms. GlobDesOpt allows selection of different design parameters, easy implementation of cost functions and models for different robotic structures. Currently, two main options are provided to optimize the design of robots in single or dualarm configuration. In this work we compare the optimization results when employing the three different solvers, along with the optimization results when optimizing the design of a surgical slave robot for the two arms independently or simultaneously.
The manuscript is thus structured as follows. Section II introduces how robot design affects the robot's model and performance, along with a brief introduction of BO, PSO, and GA. Section III presents GlobDesOpt's features and the optimization options currently implemented. Section IV shows the optimization results and, finally, conclusions are drawn in Section V.

II. PROBLEM FORMULATION
In this section, robot kinematic modelling and manipulability measure are presented, along with a brief description of Bayesian Optimization.

A. ROBOT KINEMATIC MODELLING
Different methods exist to describe the kinematic model of robots, one of the most widespread being based on Denavit-Hartenberg (DH) convention [20]. According to this method, the transformation matrix i−1 T i ∈ R 4×4 , relating the pose of a link i to the preceding one, can be described by means of 4 parameters d DH ,i , θ DH ,i , a DH ,i , α DH ,i and the joint angle q i , where d DH ,i is the distance between two consecutive joints along the current joint axis, θ DH ,i the tilting between the two joints about the current joint axis, a DH ,i the link's length, and α DH ,i the angle between the two joint axes (Figure 1). These parameters can be all stacked together generating the so-called DH table, which consists of as many rows as the number of joints n.
Since the end-effector's pose of a robot is a function of the DH parameters and of the joint values, the DH parameters and the type of joints (which are typically revolute or prismatic) highly affect the performance of the robot in executing a desired task. An optimal choice of these parameters is thus needed to guarantee the best performance.

B. ROBOT MANIPULABILTY
One typical approach to measure the capabilities of a robotic system is the dexterity. There exist various dexterity indices [21] in the literature, but they are mainly derived from the manipulability ellipsoid. An ellipsoid is defined by a core matrix H as z T (H) −1 z ≤ 1 and the singular values of the core matrix H equal the square of the semiaxes' length of the ellipsoid, whereas its determinant can be shown to be proportional to the ellipsoid volume.
For kinematic purposes, it is common to refer to the kinematic manipulability ellipsoid [22], which results by considering the set of joint velocities such thatq Tq ≤ 1. By considering that the tip twist is computed as v = Jq ∈ R 6 , with J the robot Jacobian matrix, it turns out that the kinematic manipulability ellipsoid is defined as: with H(q) = J(q)J T (q). When the robot approaches a singular configuration, the ellipsoid deforms to a line or a point (i.e. null volume ellipsoid), since at least one of the semiaxes reduces to zero. In this work we will be using the standard kinematic manipulability measure: Clearly, the robot design affects the dexterity measure by altering the Jacobian matrix.

C. BAYESIAN OPTIMIZATION
Generally, analytically defining a cost function to optimize in terms of the robot design parameters (the DH parameters VOLUME 10, 2022 and the type of joints) is intractable or even impossible. In this case, the derivatives of the function are unknown, and exploratory techniques need to be employed. Recently, in the field of machine learning, BO has found great popularity to solve optimization problems with black-box cost functions. Given a cost function f (x), BO attempts to find the global optimum in a minimum number of steps, by incorporating prior belief about f (x) and updates the prior with samples drawn from f (x) to get a posterior that better approximates f (x). Differently from other global optimization methods such as evolutionary algorithms, BO exploration is more probabilistic and exploits the model probability estimate to search through more uncertain regions. The steps performed in the optimization are: 1) Initialization: a random number of NumSeedPoints configurations are created and the cost function evaluated. 2) GP Regression: a surrogate model of f (x) is created based on the sampled points. Generally, a Gaussian Process (GP) regression model is employed. The GP outputs a probability distribution which is used to identify the sampling points. 3) Sampling: an acquisition function exploits the output of the GP to identify the sampling points, taking into account the exploration and exploitation trade-off to either evaluate points with low mean, or those with high uncertainty, respectively. The new point is added to the current dataset. Steps 2 and 3 are then repeated until the maximum number of iterations is reached. Typical acquisition functions are Probability of Improvement (PI) [23], Expected Improvement (EI) [24] and Lower Confidence Bound (LCB) [25].
In robot design, generally the numbers of parameters to optimize is not large and BO proved to be very efficient with dimensional problems with up to 20 parameters [17]. However, new approaches to deal with higher dimensional problems are being explored in the literature [26], [27].

D. GENETIC ALGORITHM
GA [28] belongs to evolutionary optimization methods, which try to find global optima of a cost function by mimicking the mechanisms of natural selection and genetic evolution. The main advantages of GA are: • it is derivative-free, thus it doesn't require to compute the derivatives of the cost function; • it finds the cost function optima in a trial and error way, by using multiple individuals per iteration; • it allows to have very large populations and, thus, effectively explore large search spaces. Due to its capabilities, GA has been widely used to identify optimal design configurations [29]. Different works employed GA for optimal robotic design [13], [30], [31].
In GA the following steps are performed: 1) Initialization: the algorithm begins by creating an initial population of PopSize individuals. This initial population can be assigned by the user or generated by means of the CreationFunction. Generally, the algorithm stops if the fitness value is lower than a specific threshold, if the maximum number of generations is reached, or if the best fitness value doesn't change for a certain number of generations (Stall Generations). At the end of the procedure, the population of the latest generation is returned.

E. PARTICLE SWARM OPTIMIZATION
Similarly to GA, PSO [32] makes no assumptions about the problem being optimized and can search very large spaces of candidate solutions, and it is derivative-free.
PSO solves the optimization by having a set of candidate solutions (the particles) and moving these particles around in the search-space according their position and velocity. Each particle's movement is influenced by its local best known position, but is also guided toward the best known positions in the search-space, which are updated as better positions are found by other particles. This is expected to move the swarm toward the best solutions. 1) Initialization: an initial swarm of particles of SwarmSize size is randomly created individuals. This initial population can be assigned by the user or generated by means of the CreationFunction. Similarly, initial random velocities are initialized too. 2) Function Evaluation: the objective function at each particle location is computed, determining the best function value and the best location. 3) Velocity and Location Update: the velocity of each particle is updated by taking into account the previous velocity, the difference between the current position and the best position the particle has seen, and the difference between the current position and the best position in the current swarm. The particles' positions are then updated with the new velocity. Steps 2 and 3 are repeated until the stopping criteria are met, e.g if the fitness value is lower than a specific threshold, if the maximum number of generations is achieved, or if the best fitness value doesn't change for a certain number of iterations.

III. GLOBDESOPT FRAMEWORK
In this section, GlobDesOpt is presented, describing its functionalities and workflow.

A. FRAMEWORK DESCRIPTION
GlobDesOpt is an optimization framework to identify the best robot design, given a specific cost function to minimize. GlobDesOpt is meant to be user-friendly and versatile, in order to be used in different scenarios and by inexperienced users. Figure 2 shows the workflow of GlobDesOpt.
Currently, GlobDesOpt allows optimizing for all the DH parameters of each robot's link d DH ,i , θ DH ,i , a DH ,i , α DH ,i and for each joint type (prismatic or revolute). The choice of what parameters to optimize for is left to the user. Additionally, GlobDesOpt supports the following design optimizations: • single arm optimization: one single kinematic chain is optimized; • dual-arm optimization: the design of two arms in a dual-arm configuration is simultaneously optimized. In this configuration, the two robots can have the same or completely different kinematic structures. In the former case, one arm is just a copy of the other. In the latter case, the DH parameters and joint types of the two robots can be optimized independently. In both cases, however, the distance between the bases of the two robots d is added as additional optimization variable. A .xml file is used as configuration file, passing information about the optimization to perform. The path to an initial DH table of the robots can be specified in the .xml file, especially if only some parameters need to be optimized. The robot model is then generated by means of a robot class function, which can also be easily edited by the user according to their individual needs. The user selects what parameters to optimize for, creating a structure of indexes for each parameter in the configuration file. Additionally, lower and upper bounds for the parameters are specified. If the dual arm configuration is chosen, with the two arms having different structures, then the indexes of the parameters for both arms need to be specified too.
Once the optimization variables and their bounds are specified, the solver is called and the optimization performed. Matlab Bayesian Optimization [33], Particle Swarm Optimization [34], and Genetic Algorithm [35] toolboxes are used in this work, and the user can specify different optimization options such as the number of iterations, the size of the population and of the swarm (for PSO and GA), the acquisition function and the exploitation/exploration ratio for BO. As described in II-C, BO runs iteratively and builds a GP model mapping the optimization variables to the cost function. By default Matlab uses ARD Matérn 5/2 kernel for the model. At each step of the optimization, the robots' structures are updated, the dexterous workspace calculated, and the cost function computed. The process is repeated until the maximum number of iterations is reached.
When the optimization finishes, the optimized parameters are returned, along with a plot of the robots' structures and workspaces.
It is worth noticing that the each solver allows including constraints in the optimization, both as bounds on the optimization variables and as additional constraint functions. These can be either added as external functions or included in the cost function itself as penalization terms. VOLUME 10, 2022 B. SINGLE ARM OPTIMIZATION Different cost functions can be implemented in GlobDesOpt, however, in this work the goal is to find the optimal design in order to improve the dexterity of a miniaturized robotic structure.
The dexterous workspace (WS) is a subset of the reachable WS, which is the set of Cartesian points that the robot can reach. In order to generate it, Monte Carlo sampling is used, randomly choosing different joint values [21]. The number of samples N s for generating the reachable workspace is userdefined.
For each configuration q n s , with n s = 1 . . . N s , the dexterity measure δ(q) defined in (2) is computed, which considers both the position and orientation terms. The user specifies an acceptance rate (α) to select only those configurations, and the corresponding Cartesian points, that have a dexterity measure above the acceptance rate. The dexterous WS is then defined as: where P is the robot's tip position, δ M , δ m are the maximum and minimum dexterity measures in the explored configurations, D q is the set of dexterous joint configurations, and D a point cloud of dexterous Cartesian points. This point cloud is then converted into a 3D volumetric shape with the Alphashape Matlab function. Once the 3D shape is created, its volume V D can be easily computed with Matlab's Volume function. As the Alphashape generates a 3D volume of triangular meshes, the Volume function adds up all the volumes of these meshes, returning the total volume of the 3D shape.
In order to take into account both the dexterity of the robot and the size of the dexterous WS for the optimal design, the global dexterity measure can be defined, similarly to [11], as: Supposing that the volume is discretized equally into N v points, then (4) can be computed as: where δ i = δ(q i ) is the dexterity value at each configuration in the workspace. The cost function for the single arm optimization can then be defined as: where the minus sign is needed as Matlab requires a function to minimize and the logarithm is just used to scale the cost.

C. DUAL-ARM OPTIMIZATION
In a dual-arm configuration, the two robots work together to a common task and it is thus important to maximize the common dexterity. To identify the common dexterous WS, first the two dexterous WS for the two arms D 1 , D 2 as in (3) are computed independently. Then, the dual arm dexterous workspace is computed as D dual = D 1 ∩ D 2 , which comprises all the points in common. This generates a new point cloud, which is then used to generate a 3D Alphashape, for which the volume V D dual can be computed. Figure 3 displays an example of the evaluation of the common dual-arm WS. Similarly to the single arm optimization, the goal is to maximize both the dexterity of the dual-arm robot and the size of the dexterous WS and, therefore, the global dual-arm dexterity measure is considered: Lee at al [10] defined the dual-arm dexterity as the volume of the manipulability ellipsoid resulting from the intersection of the two arms' manipulability ellipsoids. This measure only depends on the joint configurations of the two robots and, in their approach, it is computed locally for specific configurations in a motion task. In order to optimize the robot design and be the least task-specific, it is necessary to compute the dual arm dexterity for different possible combinations of the joint configurations of the two arms, within the common WS. Therefore, in this work, for each configuration of one of the two arms q 1,i ∈ D dual , the dual-arm dexterity measures δ dual i,j for all the configurations of the second arm q 2,j ∈ D dual can be computed. As the goal is to optimize the design in order to maximize the system's dexterity, for each configuration i in the dexterous WS the dual-arm dexterity is computed aŝ δ dual i = min j (δ dual i,j ). Finally, the global dual-arm dexterity measure considered in this scenario results to be: Algorithm 1 summarizes the computation of dual . Based on the current idea, in the case of the dual arm configuration the optimizer would lead the two robots to be for j = 1 . . . N 2 do 6: δ dual i,j = getDualArmDext(q 1,i , q 2,j ) 7: end for Compute minimum dual-arm dexterity 8:δ dual i = min(δ dual i,1 , . . . , δ dual i,N 2 ) 9:

end for
Compute average dual-arm dexterity 10:δ dual = 1 Get global dual-arm dexterity measure 11: dual =δ dual V D dual return dual 12: end function as close as possible to each other, since the common volume would then be maximized. However, the smaller the distance between the two robots, the higher the risk of collision.
In order to prevent this, a safety measure S is introduced in the cost function for the dual arm. The safety is defined as: where d is the distance vector between the two arms' bases, R max is the maximum reach, K 1 = 7, K 2 = 2 are two tunable parameters. Such function was chosen to have smooth variations in the safety values and such that min S = 0 at low distances and max S = 1 for large distances between the two arms. Note that in the dual-arm configuration, the distance vector d is added as an optimization parameter. The maximum reach is a measure to quantify how far apart the arms can be without ever colliding with each other. Once the reachable WS for each robot is obtained, each Cartesian point is projected onto the distance vector d. The maximum reach is computed as: where P 1,n s , P 2,n s (n s = 0, . . . , N s ) are the points in the reachable WS of the two arms. The minus sign is needed because the reach of the second arm is considered to be in the direction opposite to d, when all vectors are expressed in a common reference frame. Figure 4 shows an example to compute R max .
From (9), it turns out that dual arm configurations where ||d|| ≥ ||R max || are the safest, i.e. S 1, as the two arms would never collide. Smaller distances, instead, increase the chance of collisions and, consequently, lower the safety mea-  sure. Finally, the cost function for the dual-arm optimization is then defined as: where the safety measure is used as a penalization term to prevent the two arms from being too close to each other.

IV. RESULTS
In this section the design optimization results are shown, presenting a comparison of the optimal design when using the three different solvers for a dual-arm robot, along with the optimization results for a surgical slave robot both in single and dual-arm optimization mode.

A. COMPARISON OF THE OPTIMIZATIONS
In this test GlobDesOpt was employed to optimize the design of a 4-DOF dual-arm robot, with two identical arms. Table 1 reports the initial DH table of the manipulator's arms. The optimization parameters in this scenario are the first and second link lengths l 1 , l 2 , the joint types of joint 1 and 4, the distance along x, y between the two arms' base frames. The initial robot's configuration is set to be PRRR, meaning the first joint is prismatic and the others are revolute. The joint bounds are set to 0 10.0 mm for the prismatic joints and −π/4 π/4 for the revolute ones. The bounds on the link lengths l 1 , l 2 are instead set to 0 10.0 mm and 3.0 25.0 mm respectively, whereas the bounds on the distance are set to 0 15.0 mm both for x, y direction. In this exemplary test only 4 parameters and two joint types are optimized, therefore a brute force space search could be sufficient too, even though the number of computations would depend on the discretization of each parameter. In a more general case, more parameters might need to be optimized and brute force search strategies may be highly computational inefficient. Therefore global optimization approaches might be more beneficial.   All three solvers BO, PSO, and GA are here tested and their results compared. The maximum number of iterations was set to 100 and a population of 50 individuals used for GA and for PSO, whereas for BO EI was used as acquisition function and 1000 initial points generated. Additionally, the number of stalling iterations both for GA and PSO was set to 10, 50 · 10 3 points were used to generate the WS in both cases, and α = 0.6 is chosen to select the points belonging to the dexterous WS. Table 2 reports the comparison of the performances of the three solvers, and Figure 5 plots the optimized robot design and its reachable and dexterous workspaces.
In all cases, the optimal joint types result to be RRRP, namely the last joint is prismatic. BO and GA also result in same link lengths, but different base offset d. An advantage of GA and BO over PSO is that they allow specifying the optimization variables to be integers. This is of great importance because it allows taking into account constraints due to the resolution of the manufacturing process. PSO, instead, would need a proper adjustment or rounding strategy to be implemented. Given the fact that BO and GA produce the same robotic design, the volume of the dexterous WS for the two robot arms is also the same. The volume of the common dexterous WS and the average dual-arm dexterity are however slightly larger for GA optimization, due to the different relative placement of the two arms' bases. This also results in a larger safety value (0.3532 for GA and 0.2367 for BO) and smaller value of the cost function to minimize (−5.0157 and −4.7273 respectively). PSO is the solver that leads to worst performances in terms of volume of the dexterous WS and average dual-arm dexterity.
In terms of computational cost, we used a Windows x64 16 core machine with Intel i9-10980HK CPU and BO results to be the fastest solver: PSO takes 5.6 times longer and GA 6.2 times longer.

B. SURGICAL ROBOT DESIGN OPTIMIZATION
GlobDesOpt has also been employed in this work to optimize the design of a dual-arm surgical robotic instrument. Due to the satisfying and faster performances of BO optimizer from the performed simulation tests, only BO is here employed as solver. The robot consists of two identical arms, whose kinematic model is shown in Figure 6 and the DH parameters in Table 3. Each arm has 6 joints, with the first one being prismatic and the others revolute. The joint limits are set by design constraints to 0, 10.0 mm for the prismatic joint, −π, π for the first revolute joint, and −π/4, π/4 for the others. Only the three link lengths l 1 , l 3 , l 5 are here optimized, whereas the other parameters are set by design considerations.
Even though only the three link lengths are here optimized, standard optimization procedures would still be challenging to implement due to the difficulties in analytically formulating the cost function and its derivatives, especially for the dual-arm configuration.
To evaluate how to optimally design the dual-arm robotic system, two different optimization tests are performed. For both optimizations, EI was used as acquisition function, 150 · 10 3 points were used to generate the WS in both cases and α = 0.7 is chosen to select the points belonging to the dexterous WS.

1) INDEPENDENT ARM OPTIMIZATION
The design of the arms is optimized without any consideration of the dual-arm configuration. Since the two arms have an identical structure, GlobDesOpt is used with the single arm optimization option and the optimization parameters are set to be the robot's link lengths l 1 , l 3 , l 5 , with their limits set by design considerations to 3.0, 25.0 mm. In order to keep the robot small enough for surgical applications, the maximum sum of the link lengths is set to be L max = 47.0 mm and therefore a linear constraint l 1 + l 3 + l 5 ≤ L max is added to the cost function (6).

2) SIMULTANEOUS DUAL-ARM OPTIMIZATION
In this case, the interaction between the two arms is considered and GlobDesOpt is employed in the dual-arm optimization option. In this case, the same optimization parameters as in the single arm case are chosen (l 1 , l 3 , l 4 ), yet the distance between the bases of the two arms along the x direction is also added. From design limitations, the bounds for the distance are set to 0.0, 15.0 mm and, similarly to the independent arm optimization, the constraint on the maximum link lengths is added to the cost function (11).
Even though in these tests the total length constraint was added as a penalization term in the cost function, being BO an heuristic and derivative-free solver, any nonlinearity or discontinuity in the cost function would not be problematic. Figure 7 shows the optimized robot designs and Table 4 reports the optimized variables, for the independent and simultaneous optimizations. Since modern manufacturing methods ensure high precision, a resolution of 0.1 mm is allowed. Overall, the independent optimization, which does not consider any interaction between the two arms, resulted in a longer structure, with a total length of 45.6 mm versus 39.1 mm for the simultaneous dual-arm optimization.

C. OPTIMIZATION RESULTS
To compare the optimal design results of the dual-arm surgical robot, we compared three different cases by positioning the two arms obtained from the independent arm optimization VOLUME 10, 2022 FIGURE 7. The optimal designs and the surgical robot's WS from: 7) independent arm optimization; 7) simultaneous dual-arm optimization. In both cases the arms are at a distance of 14.8 mm.

TABLE 5.
Comparison of the results from the independent arm and simultaneous dual-arm optimizations for different placements d of the arms Robot 1 and Robot 2. The subscripts 1, 2 refer to Robot 1 and Robot 2 arms, V D is the volume of each arms' dexterous WS, V D dual the volume of the common dexterous WS,δ 1,2 the arms' average dexterities in their dexterous WS,δ dual the average dual dexterity in the common dexterous WS, and log(·) = log 10 (·).
at a distance of 14.8 mm (same as the optimal distance found from the simultaneous dual-arm optimization), 7.4 mm, and 1.5 mm. Figure 8 shows a comparison of the dexterous WS for the two optimizations in the different placements and Table 5 reports the values of the volumes of the dexterous WS V D 1 , V D 2 for the two arms and of the common dexterous WS V D dual , the average dexterities of the two armsδ 1 ,δ 2 and the average dual-arm dexterityδ dual , the safety value S (9) and the cost function values C (11).
Due to the two arms having the same structure, the volumes of the dexterous WS and the average dexterity measures for the two arms are the same in each test.
In all the three cases from the independent arm optimization, the volume of the dexterous WS does not change, as it is not a function of the distance between the two arms, but just of the arms' Jacobian. The common dexterous WS, instead, is larger when the two arms get closer. Having the two arms too close, however, makes the control more challenging, as the risk of self collisions increases. The proposed safety measure S prevents the arms from being too close and, indeed, it is at its minimum when the distance is set to 1.5 mm.
Because of smaller safety values, the cost function values are also larger as the distance decreases.
When comparing with the optimal design results from the simultaneous dual-arm optimization, where the distance is set to be the optimal one of 14.8 mm, the independent optimization leads to higher volumes of the arms' dexterous WS. This is also related to the fact that larger link lengths are obtained with the independent optimization.
The dexterity measures, and especially the dual-arm dexterity measure, result to be higher in the simultaneous dualarm optimization, as this approach actually considers the interaction between the two arms. Moreover, the safety measure is also larger, leading to a smaller cost.

D. SIMULATED CONTROL TASK
To further validate the best practice to optimally design a dual-arm robot, a motion control test is carried out. The two designs resulting from the independent arm and the simultaneous dual-arm optimizations, with the two arms at a distance of 14.8 mm, are compared.
The control task consists in a rough simulation of a knot tying, with the two arms following a Bernoulli lemniscate  path described by: 1 + sin 2 (t) β sin(2t) 2(1 + sin 2 (t)) 0 where P c = x c y c z c T is the center of the path, t ∈ π/2, 3π/2 for the left arm (Robot 1) and t ∈ −π/2, π/2 for the right arm (Robot 2).
In order for the path to be within the common dexterous WS and due to the different sizes of the WS resulting from the two optimized designs, the path is scaled for the two tests such that P c = 8.0 0.0 0.87 z max T mm and β = 0.4x range , where z max is the maximum height of the WS (54.0 mm for the independent arm optimization and 46.0 mm for the simultaneous dual-arm optimization) and x range is the span of the WS in the x direction (20.0 mm for the independent arm optimization and 15.0 mm for the simultaneous dualarm optimization). The two arms are required to follow the two parts of the lemniscate path keeping a constant orientation, in a synchronous and asynchronous motion (Figure 9), following the path symmetrically and antisymmetrically, respectively. Dual-arm dexterity measure values in the simulated motion tracking task for the two optimized designs. δ ave , δ min , δ max are the average, the minimum, and maximum values over the path and log(·) = log 10 (·).
Results show that the both optimization results achieve good tracking accuracy. However, the simultaneous dual-arm optimization design leads to higher values in terms of dual arm-dexterity throughout the motion (Table 6).

V. CONCLUSION
In this work we presented GlobDesOpt, a simple-to-use opensource Matlab package for optimal robot design, based on three different global optimization solvers. GlobDesOpt has here been employed for two purposes: • comparing the optimization results when using three different global optimization solvers (BO, GA, PSO); • optimally designing a surgical dual-arm robot and comparing the results between an independent arm and a simultaneous dual-arm optimization.
From the comparison of the three solvers, it turns out that GA and BO outperform PSO in terms of minimization of the cost function in the provided tests. Moreover GA and BO allow the parameters to be set as integers, thus taking into account numerical approximations due to possible manufacturing resolutions. Among the three, BO results to be the fastest approach. Additionally, results show that for a dual-arm robot, a simultaneous dual-arm optimization should be preferred, as it leads to overall higher dexterity.
Currently GlobDesOpt focuses only on articulated robots whose model is described by DH parameters. However, users can simply adapt to different structures by updating the provided robot class function, including the functions for the forward kinematics and Jacobian computation.
Future work will focus on extending the framework to other robotic structures like parallel or continuum robots, where the design complexity is higher. Moreover, an additional optimization layer will be added in order to optimally identify the best number of joints to consider in the design. We however hope that the community will benefit from this work aimed at simplifying the optimal robot design process and will contribute to use and enhance it.