SQP Optimization of 6 dof 3 x 3 UPU Parallel Robotic System for Singularity Free and Maximized Reachable Workspace

The determination of kinematic parameters for a parallel robotic system (PRS) is an important and a critical phase to maximize reachable workspace while avoiding singular configurations. Stewart Platform (SP) mechanism is one of the widely known PRS and it is used to demonstrate the proposed technique. In the related literature, GCI (Global Condition Index) and LCI (Local Condition Index) are the commonly used performance indexes which give a measure about the dexterity of a mechanism. In this work, Sequential Quadratic Programming (SQP)method is used to optimize kinematic parameters of a 6dof 3x3UPU SP in order to reachmaximumworkspace satisfying small condition numbers.The radius of mobile and base platforms and the lengths of the legs used in the platform are chosen as kinematic parameters to be optimized in amultiobjective optimization problem. Optimization is performed at different stages and the number of optimized kinematic parameters is increased at each level. In conclusion, optimizing selected kinematic parameters at once by using SQP technique presents the best results for the PRS.


Introduction
In 1947, Gough [1] established the basic principles of a mechanism with a closed-loop kinematic structure that allows positioning and orientation of a mobile platform to test tire wear and tear.In 1965, Stewart offered a very similar structure to Gough's mechanism to be used as an aircraft simulator.Ironically, the Gough platform is most often known as Stewart platform [1].In recent years, there have been many researches about PRS whose principles are originated from Gough.They have numerous advantages compared to the serial robotic systems such as high load carrying capacity and precise positioning ability.A PRS has a mobile platform connected to a fixed base by parallel legs instead of serial links.
Besides their advantages compared to serial manipulators, PRS have a small workspace and require a more complex mechanical design.Controlling the motion of a parallel mechanism is a challenging task because of its complex kinematic structure and variable dynamics through its workspace.To overcome this issue, many researchers tried to find alternative parallel structures like redundant PRS [2][3][4].In recent years, many different parallel robotic structures have been investigated for various applications.
The performance of a PRS is mainly affected by its geometrical parameters.Determining geometrical parameters is called geometric synthesis.One of the most important measures for kinematic performance is called dexterity which describes the ability of a mechanism to move and apply forces in arbitrary directions as easily as possible.However, the most suitable index to define the accuracy and rigidity performance is the condition number, which LCI (Local Condition Index) and GCI (Global Condition Index) are related to.LCI is defined as the ratio of the maximum and minimum singular values of the Jacobian matrix of a PRS [5].It varies between 1 and infinity [6,7].GCI index is a measure of dexterity and the average of LCI calculated for entire workspace.
Dexterity and workspace are the most important criteria for optimal designs of parallel manipulators.There are many researches investigating the effects of geometric parameters on dexterity, workspace, and singularity.Some of them are about maximizing well-conditioned workspace [8] using different optimization methods.Monte Carlo method [8,9], multiobjective/genetic algorithm [10][11][12], controlled random search (CRS) [13,14], analytic algorithm [15], and neurogenetic algorithm are mostly used ones [16].G. Abbasnejad offered a practical algorithm to find the singularity free cylindrical workspace of parallel manipulators using particle swarm optimization method (PSO) [17,18].In the literature, there are other performance indexes such as GSI (Global Stiffness Index) and VTF (Velocity Transformation Factor) [19].Therefore, some researchers optimized PRS for multiperformance indexes which are called MOO (Multiobjective Optimization) [20][21][22].These optimization methods are used for investigating the best geometric parameters for structures of parallel manipulators.On the other hand some researchers compared convergence performances of optimization methods [23].
In this research, we used SQP method which is Lagrangian based [24] for the optimization of kinematic parameters to reach the maximum workspace for a PRS.This paper is organized as follows.In the following section, the mechanism of a PRS and its geometric parameters and Jacobian matrix formulation for 3x3 SP are briefly described.In Section 3, our performance indexes and workspace assumption are described.Section 4 illustrates the SQP and optimization algorithm.In the last section, we discussed the results of the SQP optimization procedure.

Mechanism of the PRS
A PRS is constructed from parallel connected legs between mobile and base platforms.There are two major classification criteria of a PRS which are based on (a) joint types of connections of legs to the mobile and base platform and (b) the number of legs.The array of the types of joints are universal, prismatic, and universal, respectively.So it is called UPU.Specifically, linearly driven six legs which are connected to a mobile and a base platform cause a full spatial motion of 3-axis translational and 3-axis rotational motion.This PRS is called Stewart Platform [20,25,26] and classified by the number of leg connections on upper and lower platforms.A 3x3 SP includes 3 connections on lower and upper side of the legs.Each two of six legs share cocentered joints on each side of the legs."Figure 1" illustrates the general structure of SP that 3x3 special configuration is used for the current research.
An SP has 48 design parameters affecting its workspace and kinematic behavior [27].These parameters are spatial positions of upper and lower joints (2x18 parameters) and upper/lower leg lengths (2x6 parameters).By placing connection points on a circle, these would reduce to upper/lower circle diameters/radiuses passing through upper/lower joint centers, r p and r b .In this research, optimization parameters are these two radiuses and length of the legs,   .

Geometric Parameters of PRS.
The lengths of the legs (  ), which change according to the desired position and orientation of the mobile platform, can be calculated using the inverse kinematics equations, x (x, y, z, , , ).In this case, there are six unknowns,   (i = 1 . . ..6). "Figure 1" helps us understand the geometry of SP.Mobile and base coordinate frames, {P} and {B}, are defined.They each are linked to the mobile platform and the base.The origin of frame {P} is located at the center of a circle passing through the center of joints between the mobile platform and the legs.The Z p axis is pointing upward, and the X p axis is the bisecting line of the angle (O P P 6 -O P P 1 ).The origin of frame {B} is located at the center of a circle passing through the center of joints between the base and the legs.The Z b axis is pointing upward, and the X b axis is the bisecting line of the angle These two coordinate frames lead us to determine the position of the links attached to the base at (B 1 . .., B 6 ) and the mobile platform at (P 1 .., P 6 ). The The angles (X p , O p P i ) and (X b , O b B i ) are denoted as  i and Λ i , respectively.Therefore, we have = −    +   +    (4) "Eq." (5) "Eq." (8) describes joint positions B i of the base platform with respect to {B}.C  = cos(), S  = sin(), C  = cos(), S  = sin(), C  = cos(), S  = sin() The length of the legs can be calculated as follows: Using (4), we get length of the legs as 2.2.Jacobian Matrix.Jacobian matrix of a robotic system represents the linear transformation between the speeds in joint space, •   , and world space, ẋ , of a robotic manipulator.Jacobian matrix provides the velocity ratio from the joint space to the space of the end-effector [26].
For a parallel manipulator, the length of the legs (  ) as a vector is given as

Journal of Robotics
The position of the end-effector is x = (, , , , , )  (14) So, the relation equation between the position of the endeffector (x) and the length of the legs (  ) can be represented as By differentiating (15), we obtain Consequently J is the Jacobian matrix of the mechanism, which is used in calculating the performance indexes of the PRS, and is utilized in this research as

Performance Indexes and Workspace
3.1.Dexterity, GCI, and LCI.Dexterity is an indicator of the ability of a parallel manipulator about how rates of joints transform to the rates of the mobile platform and vice versa.It depends on the position and the structure of the parallel manipulator and can be determined by using the condition number of its Jacobian matrix as where Cond (J) is the condition number of Jacobian matrix and is defined as ‖.‖ represents the norm of the related vector or matrix.While calculating the condition number of a Jacobian matrix, different types of norms are used in the literature, and the most used ones are Frobenius and Spectral norms [1].The difference between these norms is about their calculation methods.If the Frobenius norm is considered, then we can write In this case, the Frobenius norm is defined as the extracting roots of the quadratic sum for each element of the Jacobian matrix.GCI (Global Condition Index) is defined all over the workspace as the average of Cond(J) = LCI (Local Condition Index).
If the spectral norm is introduced, the dexterity indexes can be described as where SV max(J) and SV min(J) represent the maximum and minimum singular values of the Jacobian matrix (J), respectively.The values of Cond (J), which are directly related to the singular values of the Jacobian matrix, are between one and positive infinity.All singular values will be the same in the Jacobian matrix, and the manipulator is isotropic when Cond (J) is equal to 1.
Otherwise, when Cond (J) is potentially equal to positive infinity, then the Jacobian matrix is singular.So, when the condition number is considered, the value of the condition number should be minimum for a better dexterity.
The condition number is used as an index for measuring the accuracy/dexterity of a robotic system, and also it possesses information of whether a position is close to singularity or not.However, it is not entirely possible to define a mathematical proximity to singularity for a PRS whose degree of freedom is a mix of translational and rotational motion.This is because the use of condition number is just an index and yet it has an advantage as a single number which is used to describe the overall kinematic behavior of a PRS.
Because of the complex definition of the condition number, we cannot calculate its analytical form as a function of the position parameters except for very simple robots.However, robust linear algebra tools can calculate it numerically for a given spatial position.Parallel manipulators that have both translational and rotational dof have a major drawback about condition number.Because the unit of translational motion is different from the rotational motion, the Jacobian matrix is not homogenous [28].Thus, there exist some researches about the normalization of Jacobian matrices which are partly rotational and translational.Some of them defined their own methods for normalization and some of them divided Jacobian matrix by the degree of the matrix [29][30][31][32].

Uniformity.
Uniformity is another global performance index for PRS.It is an indicator of the uniform dexterous workspace.Calculation of uniformity index is shown in equation (25).Uniformity index value changes between 1 and infinity.The index value of the best uniform dexterous workspace is equal to 1.

Workspace.
The number of possible mechanical configurations of serial robotic systems is less than that of parallel robotic system configurations.Thus, comparison of workspace of serial robotic systems is applicable.For example, one can easily assess that RRR robotic system has greater workspace capacity than a cartesian configuration because of its angular workspace capacity.However, because of the vast number of mechanical configuration possibilities of parallel robotic systems, no such comparison could be performed easily [24].One can manage to compare 6dof parallel robotic systems that are abbreviated as SP in this study.These classes of full-motion-capable spatial mechanisms are classified as 3x3, 6x3, and 6x6 based on their number of upper and lower platform connection points.So, 3x3 has triangular base and mobile platform, while 6x3 configuration has hexagonal base and triangular mobile platform configuration.Lastly, 6x6 has hexagonal base and mobile platform.Considering the first, the second, and the third configuration of similar dimensions and with the same minimal and maximal leg lengths, the third has the largest workspace, followed by the second and then the first.It must be noted that the workspace volume of the 6x6 structure is always approximately 30% greater than that of the 6x3 mechanical structure.Similarly the 6x6 configuration volume is 70% greater than that of the 3x3 workspace.Thus the 6x3 volume is approximately 25% greater than that of the 3x3 configuration workspace [1].
Additionally, the effect of geometrical parameters on parallel robotic systems is much greater than that on serial robotic systems.For instance, external stiffness over a given workspace of a 6x6 Gough Platform can differ by 700% for a change of only 10% of the platform radius,   [24].
Therefore, in order to extend the reachable workspace of a 3x3 SP mechanism, objective workspace for optimization and performance evaluation is described as follows: Therefore, total 12000 spatial, 6dof, coordinates need to be satisfied after the optimization performance.First, according to defined workspace, we calculated LCI in Frobenius norm for 12000 points in the workspace and then calculated GCI."Figure 2" represents the volume of the workspace for 3x3 SP and the inner workspace inside of it to be optimized.We used different resolution, "Δa", values for each axis of the workspace as mentioned above.For example, Δa is 0.04 m for linear motion along y-axis.

Optimization by SQP
The primary purpose of the optimization is to increase dexterity and usable workspace of the platform by changing the design parameters.Besides various optimization techniques, we used SQP method for the optimization of design parameters   ,   , and   .Min and max values of the objective design parameters are set before the optimization procedure.Another constraint is Z (LCI of workspace point).
We calculated all Z (LCI) for each point in the predefined workspace shown in "Figure 2."

Cost Function Assignments. In this research, Matlab
Optimization Toolbox (fmincon) is used to perform SQP optimization technique for the PRS.The optimization procedure is accomplished in three stages for maximum dexterity with reachable workspace: (1) One-variable optimization:   is intended to be optimized (2) Two-variable optimization:   and   are intended to be optimized (3) Three-variable optimization:   ,   ,   are optimized all together For all optimization approaches, a sub cost function is defined as follows: (i) Sub cost function for one-variable optimization is  (ii) Sub cost function for two-variable optimization is (iii) Sub cost function for three-variable optimization is (iv) Cost function for all cases is We define cost functions as described above because minimum g(x) values are investigated, which occur from the sum of minimum sub cost function, f(x), values.It means that minimum f(x) values are to be found indeed.While searching minimum f(x) values, two criteria are present.One of them is the need for minimum Z and minimum   and the second criterion is the reachable workspace.Since Z never becomes smaller than 1, the algorithm is constructed accordingly.If the targeted position in the workspace is not reachable, large value of Z is assigned not to calculate real Z value in order to eliminate the unreachable point for evaluation.Thus, unnecessary calculations finding minimum f(x) for that corresponding position are avoided.Searching for the minimum of two variables,   and   , which makes the sub cost function f(x) minimum, we multiply these variables defining the sub cost function  2 ().Thus, while minimizing the cost function, we optimize variables.We constrained Z (condition number) with a maximum value of 1000 as an indicator for dexterity.We limited leg lengths and checked if the coordinates in the workspace are reachable or not.Searching for the minimum of three variables   ,   , and   , which make the sub cost function f(x) minimum, we multiply   ,   with an average of   and (Z-1) defining the sub cost function  3 ().We optimized parameters for maximum reachable dexterous workspace.

One-, Two-, and Three-Variable Optimization Algorithm.
Figure 3 represents the optimization procedure sequence for one variable in which we have to fix other variables and define their min and max values to constrain them.Secondly, we limit the leg lengths to check if the points in the workspace are reachable or not.Later on, we define the sub cost function  1 () according to one variable or  2 () for two variables and create a limit value for selected   ().After that, we check if the workspace points are reachable or not, and we eliminate the unreachable ones.Then, cost function () is defined and condition number  for each workspace point is calculated and an optimum value for   is investigated which makes the sum of () (or ()) minimum.A similar algorithm with a little difference is performed for two-variable optimization,   and   .When optimizing for two variables, first step defined in "Figure 3" is skipped and instead of that step, we define min-max values for   as a second stage, and we use sub cost function  2 () set for two variables instead of  1 () at the fifth step.For three-variable optimization, we constrained all variables defining their minimum and maximum values.As a sub cost function, we used  3 ().

Optimization Results
3x3 SP mechanism mobile platform diameter is analysed for constant base diameter utilizing SQP technique.Mobile platform/base platform diameter optimization for better workspace reach and dexterity is performed.Condition number variation for different mobile platform diameters in selected range, minimum of condition number for base/mobile platform diameter change, and the number of reachable points for different base/mobile platform diameter combinations are discussed.
In Figure 4, horizontal axis refers to scanned   values, and vertical axis refers to minimum condition number which we get from calculations without homogenization of Jacobian  matrix.We get low values of  min while  p is increasing.We are searching optimum values for most efficient workspace volume which has reachable points with low condition numbers.As a result, high mobile platform diameter refers to calculated minimum condition numbers.
In Figure 5, horizontal axis refers to scanned   values.Left vertical axis refers to scanned   values and right vertical axis refers to the mapping of condition numbers which we calculated in Frobenius norm for predefined reachable workspace.As seen in Figure 5, minimum condition numbers occur with high values of   and   which mean dexterity is getting better while they are increasing.However, it is obvious from the plot that comparable decrease exists of the mobile platform diameter with respect to base diameter to obtain minimal condition numbers.Also, the rate of this decrease is getting smaller when the base diameter is a lot higher than mobile platform diameter.The horizontal tendency of platform diameter is apparent for higher base diameters.
In Figure 6, horizontal axis refers to scanned   values.Left vertical axis refers to different   values, and the right vertical axis refers to the mapping of reachable points in the workspace which have condition numbers limited with an . Low mobile platform diameters with high base diameters seemed to result comparably in small cost function values, which means better dexterity in motion.In Figure 7, minimum cost function value is 4.8 x 10 7 as indicated in bold areas.In those areas, cost function values are minimum with optimized   and   values as shown and are compatible with Table 2.While optimized   has minimum values, cost function values are minimum.As mentioned previously, we are searching for minimum cost function value for better dexterity.In Table 2, there is not one optimum solution for the minimum cost function value.For example,  , ,  , = 0.074 [m], 0.129 [m] has minimum cost function value of (4.8xE7) and also  , ,  , = 0.077 [m], 0.125 [m] has minimum cost function value of (4.8xE7).
For better accuracy, we increase the resolution in solution area in Table 3.For the purpose of interpreting the results, we used cost function values.However, the condition number Optimum values are in Table 3. Minimum cost function value is 3.2xE7.Solution set with minimum cost function   4, our constraints are as in Table 5.For three-parameter optimization, we used  3 () as a sub cost function and obtained results in Tables 4 and 6 using the optimization algorithm as shown in Figure 8.
For Table 7, we used optimized values from Table 3, and we generated Table 7 using inverse kinematic equations.
Table 7 presents L max , L min ,  (uniformity index), and reachable points for optimum radius values of mobile and base platforms.According to these calculated values using optimized parameters from Table 7, GCI (Global Condition Index) and LCI (Local Condition Index) values are computed for each solution.Figure 9 represents LCI values for the first  In detail, first solution is the best for the amount of the reachable points.Second one is the best solution according to GCI.And the third one is the best for the uniformity.

Conclusion
Determining the values of geometrical parameters is very important for efficiency of a PRS design.In literature, there are many researches for the optimization of kinematics of a PRS.In these works, various optimization methods have been used.In this study, we optimized 3x3 SP mechanism using SQP optimization method for predefined workspace and checked the efficiency of the optimization for three parameters which are diameter of mobile and base platform together with lengths of the legs.It is shown that the optimization of multi-parameters is not easy to formulate concerning cost functions and results must be evaluated by calculating the dexterity indexes and the amount of reachable points.Firstly, we optimized   and then, we optimized   and   simultaneously.Lastly,   ,   , and   are optimized at the same time.During optimization, tolerances, procedure steps, and applied/defined cost functions affect the convergence efficiency.Results show that initial values for each parameter at the beginning of the optimization determine the tolerance of the convergence for the optimization.As a further research, various optimization methods can be performed while increasing number of kinematic parameters optimized searching different cost functions.  : Radius of a circle passing through center of joints between base platform and legs   : Radius of a circle passing through center of joints between mobile platform and legs   : Angle between connection points of legs on base platform   : Angle between connection points of legs on mobile platform R: Rotational transformation matrix.
angle (O p P 2 , O p P 3 ) is denoted as   .Angles (O p P 1 , O p P 3 ), (O p P 3 , O p P 5 ), and (O p P 5 , O p P 1 ) are each equal to 120 ∘ .Similarly, angle (O b B 1 , O b B 2 ) is denoted by   .Angles (O b B 1 , O b B 3 ), (O b B 3 , O b B 5 ), and (O b B 5 , O b B 1 ) are, also, equal to 120 ∘ .Now, the length of the legs (  ) can be calculated for  = 1, 2, . . ., 6

Figure 2 :
Figure 2: Volume representation of the workspace to be optimized for the SP and the overall workspace.
(i) P x = Changes from -0.06 to +0.06 m in 0.04 m regular intervals (ii) P y = Changes from -0.06 to +0.06 m in 0.04 m regular intervals (iii) P z = Changes from +0.3 to +0.4 m in 0.02 m regular intervals (iv) , ,  = Changes from -5 to +5 in 2.5 ∘ regular intervals

Figure 4 :
Figure 4:   (min condition number) variations for different   values while   is fixed to 0.175 [m].

Figure 5 :
Figure 5:   with respect to   and   in predefined reachable workspace.

0
.050 <   < 0.180 0.070 <   < 0.180 0.080 <   < 0.180  , = 0.075    = 0.075    = 0.088 0.090 <   < 0.180 0.100 <   < 0.180 0.110 <   < 0.180    = 0.090    = 0.100    = 0.110 0.120 <   < 0.180 0.090 <   < 0.250 0.050 <   < 0.250    = 0.120    = 0.090    = 0.075 upper bound of 1000.As indicated in Figure 6, while   values decrease,   values increase, and numbers of reachable points elevate.The upper bound of Z value is equal to 1000.Thus, the numbers of reachable points are high for low mobile platform diameters and high base diameters.During optimization process of mobile platform diameter, different initial values of   with various constraints are used, as shown in Table 1.While   = 0.175 [m], minimum limit of   is kept 0.050 [m] and leg lengths are scanned between 0.300 [m] and 0.450 [m].As a result, the optimum value of   is calculated to be 0.075 [m].When minimum limit for   is 0.080 [m], optimum value is found as   = 0.088 [m].As mentioned above, different initial values of   are used for optimization.Consequently, when the constraint value of   changes, the optimum value of   changes.This is because, in different workspace regions, we have different local minimums.If the optimization procedure in Figure 8 is analyzed thoroughly, while the lowest constraint of   is 0.050 and 0.070 [m], we selected optimization scanned range between 0.050 and 0.180 [m].Then, we obtain an optimum value of   0.075 [m].When the min constraint value of   is elevated to 0.080 [m], the first local minimum of   0.075 [m] is passed.Optimized value of   shifts from 0.075 [m] to 0.088 [m].Out of nine different optimization combinations, one can understand that we have two local minimums which are 0.075 [m] and 0.088 [m].Then, combinations that include the initial value of 0.090 [m] and larger will result in an optimized value of start point of the scanned range, 0.090 [m], 0.100 [m], 0.110 [m], etc.Consequently, we found two local minimum values (0.075 [m] and 0.088 [m]) for   having maximum reachable points which are limited to upper bound condition number,
, and 11 are calculated from the homogenous Jacobian matrix and these calculated values are divided by the dimension of Jacobian matrix for homogenization.

Table 1 :
Optimization of mobile platform diameter for various range constraints.
Z < 1000 Amount of Reachable Points in Workspace

Table 2 :
Optimizing mobile platform diameter,   .Z, of 1000.Optimization results for two variables   and   and their cost function values with different initial values of the optimization procedure are as shown in Table 2.For example, while the initial value of   is equal to 0.125 [m] and an initial value of   is equal to 0.070 [m], cost function value is equal to 4.8xE7, optimum value of   is equal to 0.074 [m], and optimum value of   is equal to 0.12 [m].Optimum values of diameters   and   depend on constraints.While generating Table 2, constraints are r b : 0.125 [m] < r b < 0.175 [m] and constraints for

Table 3 :
Optimization of mobile platform diameter,   , by narrowing constraint target area.

Table 5 :
Constraints for three-variable optimization of the PRS.