Root Based Optimization Algorithm for Task-Oriented Structural Design of a Multi-Axial Road Test Rig

A novel root based optimization algorithm (ROA), which mimics the proliferation of plant roots, is proposed to evaluate the performances of a multi-axial road test rig over a prescribed workspace. Plant roots have developed complex behavior patterns to search for nutrients in the soil and even exhibit swarm intelligence. The growing roots can adapt their actions such as elongation, branching, and tropistic movement according to the environment. During the process of foraging for water and nutrients, the adaptive growth behavior of roots can be simulated to tackle the optimization problem. The efficiency of the proposed ROA has been validated by comparing it with well-known intelligent algorithms using classical and IEEE CEC 2019 benchmark functions. Computational results indicate that ROA outperforms other algorithms to search a global optimum on several benchmarks. Furthermore, the results are tested with nonparametric statistical methods, e.g., Wilcoxon rank-sum test. Besides, ROA is applied to the optimal dimensioning synthesis of the multi-axial road test rig (MRTR) based on the orthogonal 6-RSS parallel manipulator to promote the device’s performance. The inverse kinematics of the multi-axial road test rig based on the asymmetrical orthogonal 6-RSS parallel mechanism is analyzed in detail. Applying the achieved dimensionless Jacobian, the performance criteria involving the global isotropy index, the force/velocity transmissibility index, the decoupling index are established. By simultaneously taking account the proposed criteria of performance, structural optimization is carried out resorting to ROA. Finally, a set of optimized dimensional parameters is obtained. The practical application further illustrates that ROA has the potential to deal with complex optimization problems.


I. INTRODUCTION
Many nature-inspired intelligence algorithms have been developed to address complex and nonlinear problems, which cannot be solved effectively by traditional methods. The existing meta-heuristic algorithms may be divided into two main categories: Evolutionary algorithms (EA) and Swarm intelligence algorithms (SI) [1].
Evolution is a stochastic process and emerges at the population level. A generation of individuals will bring forth offspring by recombination and mutation, and the fittest are selected to survive each generation [1]. Notable evolutionary The associate editor coordinating the review of this manuscript and approving it for publication was Huaqing Li . algorithms include the Genetic Algorithm (GA) [2], Differential Evolution (DE) [3], Invasive Weed Optimization (IWO) [4] and Improved Genetic Immune Algorithm (IGIA) [5], etc.
Swarm intelligence is the collective behavior of the biological population and the local interactions between individuals lead to global behavior. Most species show swarm intelligence and the well-known examples consist of bird flocks, fish schools, herds of quadrupeds, the colony of social insects and bacteria molds [1]. Therefore, inspired by these diverse swarm systems, a plethora of algorithms, such as particle swarm optimization (PSO) [6], PSO with inertial weights (GPSO) [7], Cuckoo Optimization Algorithm (COA) [8], Harry's Hawk Optimization Algorithm (HHO) [9], Grey Wolf Optimizer (GWO) [10], Whale Optimization VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Algorithm (WOA) [11], Squirrel Search Algorithm (SSA) [12], Bee Colony Optimization (BCO) [13], Fitness Dependent Optimizer (FDO) [14], Moth Flame Optimization Algorithm (MFO) [15], Slime Mould Algorithm (SMA) [16], etc, have been created. An elaborate literature review on nature-inspired algorithms is presented in [17]. In contrast to animals, plants based algorithms are scarce and only a few studies inspired by plant reproduction and root growth have been carried out. Mehrabian and Lucas [4] introduced the evolution-based Invasive Weed Optimization (IWO) through imitating the behavior of weeds in occupying territory and colonizing. The IWO is further enhanced by combining chaos theory [18]. Plant propagation algorithm (PPA) is based on the way the strawberry plant propagates using runners, which can explore the search space with long runners while short runners enable it to exploit the local space [19]. Similarly, Runner Root Algorithm (RRA) mimics the behavior pattern of strawberry and spider plants which propagate through their runners for global search and produce roots to search for local resources simultaneously [20]. Flower Pollination Algorithm (FPA) is proposed to simulate the features of flower pollination, combining local and global searches through self-pollination and long-distance cross-pollination [21]. Inspired by local and global seed dispersal methods of trees, Forest Optimization Algorithms (FOA) has been proposed [22]. Similarly, Natural Forest Regeneration Algorithm (NFR) searches the global and local space through random seed dispersal and decreasing the radius of seed dispersal [23]. Artificial plant optimization algorithm (APOA) is a methodology mapping the growing process of a tree into optimization, including three phenomena: photosynthesis operator, phototropism operator and apical dominance operator [24]. Also, Cheraghalipour [25] developed the Trees Growth Algorithm (TGA) inspired by trees competition for acquiring light and foods. According to the self-defense mechanisms of plants, an optimization algorithm is presented based on the predator-prey model [26]. Apart from propagation, plants have evolved complex root systems to forage for water and nutrient. Root Mass Optimization (RMO) has been developed based on root re-growing and root branching [27]. Furthermore, Artificial root foraging optimization (ARFO) is proposed which uses root branching and elongation for local and global search, and random growth of lateral roots keeps the diversity [28].
When a root is located on a fertile circumstance, it will grow faster and produce more collateral roots, otherwise, it will generate fewer branches or even cease growing. According to the No-Free-Lunch (NFL) theorem [29], there is no algorithm for solving all optimization problems. Therefore, it is still a challenge to design better optimizers. Inspired by the root proliferation strategies, a new root optimization algorithm (ROA) is developed for global optimization. Unlike RMO and ARFO, the proposed ROA adopts three new strategies: (1) emission of primary roots, (2) branching, (3) re-growth. In the emission phase, let the main roots elongate to move toward a better food source. According to the fitness, the main roots produce new roots to exploit around the local area. In the re-growing phase, let some of the branch roots randomly move to reach new places. Depending on the fitness value, the new main roots are selected from all populations. The performance and efficiency of the ROA have been tested on some benchmarks. Besides, ROA is applied to the optimal design of road test rig based on asymmetrical orthogonal 6-RSS parallel manipulator with a specific task. The kinematic and static instantaneous performance indices based on the scaled Jacobian matrix have been specifically developed for the service load simulation in this article, containing force isotropy index, force/motion transmission index and decoupling index.
The major contributions of this article are as follows: 1. A novel plant roots-inspired intelligent algorithm (ROA) is proposed. The growth mechanisms of roots are studied thoroughly and modeled mathematically including emission of primary roots, branching and re-growth, which combines swarm intelligence and evolutionary strategies to trade-off in both exploration and exploitation abilities.
2. The proposed algorithm is validated on 29 benchmark functions including the classical as well as modern CEC 2019 benchmarks. The efficiency is also studied by rigorous comparison with existing well-known nature-inspired optimization algorithms through nonparametric statistical analysis.
3. The inverse kinematics of the multi-axial road test rig based on the asymmetrical orthogonal 6-RSS parallel mechanism is analyzed in detail and the Newton-Raphson algorithm has been proven to be efficient for solving the problems of the inverse kinematics. Based on the four performance indices, the effectiveness of ROA is investigated for the optimal dimensioning synthesis of the multi-axial road test rig.
The remainder of the paper is organized as follows. Section II describes the principle of root optimization algorithm (ROA) and validation and analysis have been performed. Section III analyzes the inverse kinematics of a multi-axial road test rig. Section IV establishes the kinematic and static performance indices. Section V examines the effect of each design parameter on the performance indices and presents the dimensional optimization with the ROA algorithm. Finally, Section VI presents the conclusions of the paper.

II. ROOT OPTIMIZATION ALGORITHM A. MATHEMATICAL PRINCIPLE OF ROA
Foraging for limited resources, plant roots develop complex behavior patterns and even exhibit swarm intelligence like animals through fast electrical and chemical signals interaction. The roots can sense the environment to adapt their actions, e.g., the growth direction adjustment is based on the mutual distance and roots can be all directed towards areas with optimal water and nutrient gradients in the soil [30]- [32].
Several researchers have investigated many models in accordance with experimental data to mimic growing-roots patterns. The roots can be divided into two subsystems: main roots and collateral roots. The key principles are composed of three phases: (1) emission of primary roots, (2) branching, (3) re-growth. The enormous random proliferation of the main roots in areas free of other roots will dominate the initial stage, providing a highly explorative capability. As iterations surge, main roots slow down their growth rate and more collateral roots manage to navigate toward the source of the resource by local fluctuations. Furthermore, current resources can be exhausted rapidly and new research must be undertaken.
Assume the plant produces M main roots, i.e. X = (X 1 , X 2 , · · · , X M ), each of which searches in the Ddimensional growth space. The i − th root is described by where the D is the number of design variables, and LB j and UB j are the lower and upper bound of growth space for j − th design variable x ij , respectively. The element of the vector is initialized as follows where rand() is a random number distributed uniformly in [0,1]. The roots in each generation are sorted in descending order (or ascending order) according to the fitness asX = (X 1 ,X 2 , · · · ,X M ), of whichX 1 andX M are the best position and the worst position, respectively.

1) EMISSION OF PRIMARY ROOTS
At first, the branch roots growth mechanism is deduced as follows: The grow radius of the i − th root is defined as where c 1 is a random vector distributed evenly in [0,1]. When the main roots reside in better conditions, it will produce more collateral roots. The number of branch roots as offspring generation is a linearly decreasing function concerning the fitness value of the main roots. Additionally, the number b of branch roots will diminish dynamically with the shrink of the grow radius using a nonlinear function.
The k − th branch root produced by i − th main root is written aŝ where c 2 is a random vector distributed evenly in [−1,1], and b = b max −i+1, where the maximum number b max of branch roots satisfies where b max 0 is the initial maximum number of branch roots, m u0 is proliferating factor and Iter and Iter max are the iteration number and its maximum value, respectively. Inspired by [33], the m u is defined as the probability of generating branch roots by main roots, as shown in Figure 2.

2) BRANCHING
The grow radius is narrowed over time and the algorithm will converge quickly to the local optimum solution. To explore the growth space sufficiently, roots grow freely similar to the method in [27] as followŝ where c 3 is a random number distributed evenly in [−1,1] and γ is defined as where c 4 is a random vector distributed evenly in [−1,1].

3) REGROWTH
After the current resources have been exhausted, new research will start to search for a new area. To enhance VOLUME 8, 2020 exploration capability, the branches are relocated through the following equation: where lévy distribution can provide more efficient exploration. The lévy flight is stated as [12] Lévy(β) = 0.01 × r a × σ where r a and r b are random numbers normally distributed in [0,1], β is a constant (β = 0.5), and σ is calculated as where (x) = (x + 1)!.
After the update of the position of the roots, the new main roots are reselected from the merged set according to the survival of the fittest and the number of new main roots in the next generation remains unchanged. The pseudo-code of the ROA algorithm is illustrated in TABLE 1.

B. BENCHMARK TEST
The ROA algorithm has been implemented in MATLAB TM computer program and its efficiency was validated with some well-known benchmark functions [34] as shown in TABLE 10. TABLE 10 summarizes the three types of test functions, the number of variables, the range of variation, and the optimal values. ROA was compared with GPSO, GA and the other state-of-the-art swarm-based IWO, GWO, MFO, and WOA algorithms.
The selection reasons for these algorithms are: 1) They are based on different search strategies involving EA and SI, and inspired by diverse species. 2) All of them are well cited and more mature. 3) They have been applied to tackle real-world problems. 4) These algorithm implementations can be obtained publicly. All algorithms parameter settings remain the same as ones in their original papers, as shown in TABLE 2.
For all the algorithms, population size and maximum iteration equal to 30 and 500 have been utilized. For each benchmark function, the ROA algorithm was run 30 times the same as the other algorithms.

1) CLASSICAL BENCHMARK TEST FUNCTIONS VALIDATION
Test functions involve three sets: unimodal, multimodal and composite functions, which are used for testing exploitation level, exploration level and practical problem-solving ability of the algorithm. Statistical results (i.e. average cost function and corresponding standard deviation) are presented in TABLE 3.
The unimodal functions TF1-TF7 have only one global optimum and can be used to evaluate the exploitation capability of the algorithm. It can be observed from TABLE 3 that ROA is very competitive with other algorithms. Obviously, ROA outperforms other all algorithms for all functions TF1-TF7 except TF5. Therefore, the proposed algorithm can provide very strong exploitation.
Unlike unimodal functions, multimodal functions TF8-TF13 are composed of an exponential number of local optima concerning the design variables size. Hence, this kind of test functions is more suitable to assess the exploration level. The results reported in TABLE 3 indicate that ROA shows better results for functions TF8-TF11 and provides suboptimal results for TF12-TF13. The results confirm that ROA can explore the growth space extensively and is capable of superior exploration capability. Composite benchmark functions TF14-TF19 with abundant local optima can test the tradeoff ability to handle complications existing in real-world problems. As shown in TABLE 3, ROA provides better results than other algorithms except for TF18 and still had a relatively comparative result in TF18.
As a summary, the results demonstrate high exploration and exploitation characteristics of the proposed ROA algorithm.

2) CEC-C06 2019 BENCHMARK TEST FUNCTIONS VALIDATION
The CEC 2019 benchmark test functions [35] in TABLE 11 for a single objective optimization problem are also utilized to evaluate ROA. As shown in TABLE 4, ROA outperforms other algorithms except in CEC02, CEC03 and CEC06 benchmarks. Even though GPSO has better results for CEC02 and CEC03 benchmarks, ROA has only a small discrepancy.

3) STATISTICAL TESTS
Three nonparametric statistical tests involving Wilcoxon test, Kruskal Wallis test, and Friedman test [36] are selected to analyze the performance of these algorithms.

III. INVERSE KINEMATICS MODELLING OF THE MRTR
The parallel manipulators possessing several advantages such as higher accuracy, larger load capacity over serial robots [37], [38] have been applied widely in various industrial fields. As a typical example in the automotive industry, the Multi-axial Road Test Rig (MRTR) is used to replicate the real road condition in the laboratory and has become an imperative instrument for vehicle durability testing [39]. The increasingly stringent requirements of service load simulation for vehicles put high demands on the performance of these testing devices. Since performances of the parallel mechanism are very sensitive to their dimensioning [40], VOLUME 8, 2020   the optimum structure parameter design plays a key role in the performance improvement of the MRTR.

A. STRUCTURE DESCRIPTION
The concrete structure and schematic diagram of the MRTR are shown in Figure 3 and Figure 4, respectively. In essence, the MRTR is a specific 6-RSS parallel manipulator as shown in Figure 3, composed of an end-effector, six coupling links, six hydraulic actuators, and seven bell cranks. Driven by six hydraulic actuators, the end-effector can exert vertical, lateral, and longitudinal forces as well as camber, steer, and brake moments at each spindle of the vehicle, as shown in Figure 4. The light-weighted coupling link and bell crank can lead to a reduction of the inertia of the moving links. Besides, the actuators are allowed to be located close to the base which also reduces the inertia of the moving rods.
The MRTR with unique topology possesses superior properties, e.g. the orthogonal configuration of limbs reduces the cross-couplings between sub-chains and complexity of the control system [41].
As shown in Figure 4, all parts connect each other with S (spherical) joint, R (revolute) joint and actuated P (prismatic) joint by a hydraulic actuator. Thus, the end effector is attached to the fixed base through six limbs, where the limb1, limbs 4-6 have the identical structure. Limb 2 and limb 3 are merged and form a 2-DOF hybrid kinematic chain, which can execute vertical or brake input independently. The brake rotation is transmitted from the bell crank to the mobile platform by the parallelogram two-bar-linkage. The brake event occurs occasionally during the vehicle durability testing, so the 2-DOF serial sub-chain plays a significant role. And its principle is illustrated in Figure 7 in detail.

B. COORDINATES DESCRIPTION
To facilitate kinematic analysis, a fixed Cartesian frame O b − x b y b z b and a moving frame O p − x p y p z p are assigned on the fixed base and at the centroid of the moving platform in Figure 5, respectively. The local reference frames a i − x i y i z i (i = 1, 2, · · · , 6) are attached to the limb i. Define   of Euler angles can be expressed as [42] where s (·) = sin (·) , c (·) = cos (·) and the same in the following. As shown in Figure 5, in the reference frame {b}, the position of the point p i on the platform can be described as  where R is the transformation matrix from the moving frame {p} to the fixed frame {b}, following Z − Y − X convention. The matrix R can be written as (15), as shown at the bottom of the next page.

C. INVERSE DISPLACEMENT SOLUTION
As shown in Figure 6 and Figure 7, the constrained equation for ith(i = 1, · · · , 6) limb can be written as where l i represents the vector of the ith coupling link in the fixed frame {b}, a i is the position vector of the ith revolute joint a i and b 1i , b 2i denote the position vector of the ith bell crank transmission arm and actuation arm, respectively. Define unit position vector of the ith transmission arm and actuation arm as s 1i and s 2i respectively, so the b 1i and b 2i can be rewritten as VOLUME 8, 2020 where b 1i and b 2i represent the length of transmission arm and actuation arm, respectively and the unit vector s 1i and s 2i can be derived as where θ i and β i are the rotation angle and constant angle of the bell crank, respectively. The Eq. (16)-(18) depicts the holistic kinematics of MRTR and reformulate them in matrix form as follows: where 6 ] T and s = [s 11 · · · s 16 ] T are the matrix of position vectors about coupling link, the joints on moving platform, joint a i and bell crank transmission arm, respectively. And S represents the coupling coefficient matrix, which is According to the Eq. (23), the function concerning θ i can be derived as: where θ = [θ 1 · · · θ 6 ] T is the matrix of bell crank rotation angles. Define s 1i = ds 1i dθ i as derivative concerning θ i and s = [s 11 · · · s 16 ] T is the matrix of the derivative function.
Based on the Eq. (25), the matrix form can be obtained where • means the Hadamard product (i.e. the elements multiplication). By resorting to the Newton-Raphson algorithm, θ i can be solved on-line with the Eq. (27)- (30) and Eq. (28) can be solved by QR decomposition to obtain the θ k .

D. VELOCITY ANALYSIS
The velocity of the joint point p i (i = 1,2, · · · ,6) by taking the derivative of Eq. (14) can be obtained aṡ (34) whereṗ b i is the velocity of p i and I ∈ R 3×3 is a unit matrix, and J piq represents a Jacobian matrix between the platform and joint point velocity, J piq = I −Rp p i × I . Differentiating Eq. (23) concerning time, the velocity function can be given by whereθ = [θ 1 · · ·θ 6 ] T = [ω b1 · · · ω b6 ] T represents the scalar angular velocity vector of the bell crank.
where J ωbq denotes the Jacobian matrix between scalar angular velocity and platform velocity. And then, as shown in Figure 8, the angular velocity ω bi of the ith bell crank can be calculated as where J ωbiq denotes the Jacobian matrix between the ith bell crank angular velocity and platform velocity, and e represents the unit vector, which is where e x and e y are the unit vector along with the x b axe and the y b axe of the fixed frame {b}, respectively.
Assume v di be the S joint point velocity of the ith bell crank actuation arm, which can be described as where· denotes the skew symmetry matrix, J vdiq denotes the Jacobian matrix between the S joint point velocity of the ith bell crank actuation arm and platform velocity.
Let v dni be the ith actuator scalar sliding velocity, which can be derived as where v dn = [v dn1 · · · v dn6 ] T denotes the sliding velocity of the actuator and J = J vdn1q · · · J vdn6q T is the kinematic Jacobian matrix mapping Cartesian velocity to joint rate.

IV. KINETOSTATIC PERFORMANCE INDICES FORMALIZATION
Performance synthesis for a multi-DOF mechanism is an application-dependent task [43] and the desired performance of MRTR consists of good kinematic conditioning, high forces transmissibility, and low motion cross-coupling. Performance evaluation indices are a prerequisite for the dimensional synthesis. The Jacobian matrix contains significant information concerning the kinematic and static characteristics of the mechanism so a variety of indices are derived from it to evaluate the kinetostatic performance of a parallel mechanism [40]. Isotropy indicates the directional dexterity degree of a mechanism and is one of the common performance measures [44]. Salisbury and Gosselin used the condition number of Jacobian to define the Local Conditioning Index (LCI) [45] and the global conditioning index (GCI) [46] for mechanism isotropy optimization.
The homogeneous Jacobian can also be accomplished by multiplying scaling matrices [47]. Applying the task space and joint space scaling matrices, Stocco et al. [48], [49] designed the global isotropy index (GII) to obtain a position independent worst-case optimum over a workspace.
As an alternative to condition number, Olds [50] proposed p, q-norm indices to quantify local and global worst-case kinematic and force performance. Lou et al. [51] applied the minimum singular value of Jacobian in the prescribed workspace to evaluate different performances. Shao et al. [52] adopted the matrix orthogonal degree of Jacobian to evaluate the motion/force transmission performance. Apart from the Jacobian-based indices, motion/force transmissibility performance indices based on the transmission angle [53] or the screw theory [54] have also been developed and are especially prevalent in the analysis of lower-mobility parallel mechanism. When a parallel manipulator possesses kinematics decoupling characteristics, it can simplify the kinematics analysis, and enable calibration and control easier, especially for the high-frequency vibration application [38]. The parallel road test rig requires minimal joint cross-coupling with good force manipulability [55]. Jin and Yang [56] classified the kinematic decoupling into the strong coupling, complete decoupling and partial decoupling. Jin et al. [57] introduced the concept of group decoupling. Chu and Gao [58] proposed the distribution difference of nonzero elements in Jacobian matrix to measure the kinematic coupling degree of forging manipulator. Shen et al. [59] employed a systematic method to analyze the parallel kinematic structures with a lower coupling degree. Legnani et al. [41] developed some new 6-DoF parallel mechanisms decoupled in one specific configuration by modifying the leg structure. Lee et al. [43], [60] put forward the local kinematic cross-coupling index (LKCI) with the angle of interaction between column vectors of the Jacobian to signify joint dependency. In [42], kinematics and VOLUME 8, 2020 statics cross-coupling indices were also explored in detail with the same mathematical principle.
Establishing the kinetostatic indices allows us to quantify the performance of the test rig. The multi-criteria based design can provide sufficient control on the range of the design parameters so it is necessary to establish multiple performance indices [61]. For the multi-axial road test rig, the dexterity, velocity/force transmissibility and decoupling degree are more concerned. Therefore, four corresponding indices are given to measure these performances in the following part.
When a specific desired performance goal exists, the Jacobian matrix can be normalized with scaling factors to remove the inhomogeneity limit [49]. Because the multi-axial road test rig has the unambiguous task requirements as shown in TABLE 7, the task-specific design approach is more suitable and practical. As a consequence, four performance measures are all developed based on the scaled Jacobian.

A. FORCE MANIPULABILITY MEASURE
Assuming the gravitational forces, inertia forces, and frictional forces unconsidered, the mapping between the generalized force F exerted on the end-effector and the joint force f is (42) According to the maximum desired task-space forces and the maximum force capabilities of actuators, the scaling matrices are given by The normalization of force transformation Jacobian J F can be derived as follows where F and f are the percentage vectors of maximum task-space forces and maximum actuation forces, respectively. The workspace-inclusive global isotropy index is a more stringent measure than the condition number, which is defined as [62] η h = GII (ξ ) = min It can be seen that a larger GII value promises a better control dexterity in general.
For the multi-axial road test rig, exerting larger generalized forces on the spindle of the vehicle is expected for the given actuation forces. Here, the minimum singular value σ min J T F is used to characterize the force manipulability of the mechanism, then the force transmission index κ f can be expressed as In a real implementation, usually the maximum task-space force is expected to be produced in all points in the prescribed workspace for the test rig. The global force index (GFI) of a defined workspace W is defined as A larger GFI value allows a better force transmission rate.

B. MOTION MOBILITY MEASURE
For succinctness, let v = v dn and the mapping between task domain velocity and joint velocity vector can be reformulated asq Similarly, the corresponding scaling matrices for velocities are The scaled velocity Jacobian J v can be derived as where q and v are the percentage vectors of maximum task-space velocities and maximum actuation velocities, respectively. Still, based on the minimal singular value, the velocity transmission index κ v can be mathematically expressed as In the automotive testing industry, the velocity of the multiaxial test rig should reach a desired value in the prescribed workspace and high-speed are always desirable. The global velocity index(GVI) over a defined workspace W is defined as A larger GVI value indicates a better velocity transmission rate.

C. CROSS-COUPLING MEASURE
Kinematics decoupling means that one motion of the endeffector only corresponds to the input of one limb or one group of limbs [56]. It is difficult to design complete decoupling mechanisms for parallel manipulators' overall workspace, but it is possible to reduce the coupling in one particular configuration or the limited workspace through optimization design.
The cross-coupling between joint motions can be interpreted as the degree of orthogonality between column vectors of the Jacobian matrix [43]. A quantitative measure of kinematic cross-coupling between two sub-chains is defined as [60] κ ij = sin α ij (61) where α ij denotes the angle between ith and jth column vectors of the Jacobian, which is Since the angle of intersection between two vectors is in the range of (0, π), it can be deduced that 0 ≤ κ ij ≤ 1.
At one specific configuration, the more general local kinematic cross-coupling index for analyzing all joints dependency is defined as [42] When κ cp = 1, each two column vectors of the Jacobian matrix are mutually orthogonal and all joint motions are completely decoupled. The global kinematic cross-coupling index (GCPI) in the prescribed workspace W is given by Since there exists no analytical means of calculating the Jacobian, numerical integration is required to determine the index value [38]. Thus above integral of global performance indices may be approximated by a discrete sum where m is the integration points of workspace. To reduce the fabrication cost and simplify the design procedure, some constraints are put forward as follows

B. APPLICATION TO KINEMATIC PERFORMANCE EVALUATION
To make the number of independent variables more manageable, the effect of individual candidate design parameters on the performance indices should be clarified.
An initial value of each design parameter is assigned as shown in TABLE 8 and a large candidate range of variation VOLUME 8, 2020 The change rule of global performance indices is shown in Figures 9-11. In Figure 9, it is clear that increasing the coupling link length leads to an accompanying increase in all global performance indices. When scale-out the initial size, change of indices inclines to be stable, indicating that the effect of coupling link length on the performance of mechanism can be ignored after the coupling link size reaching a certain length. The variables l 2 and l 4 have more effects on GFI and GVI. And, GCPI changes to a similar extent with the variation of variables l 1 , l 2 and l 4 . In conclusion, the performance will improve with the increase of coupling links length.
In Figure 10, the larger the size of the bell crank transmission arm is, the better the isotropy, velocity transmission and   decoupling become. While the GFI becomes smaller with the increase of variable b 11 , i.e. force transmission performance becomes worse. It also can be seen that with the increase of variable b 11 , four indices tend towards steady which is consistent with the change in length of the coupling link. When the length of the bell crank transmission arm reaches a VOLUME 8, 2020 certain extent, there will be no longer a significant influence on the performance of the test rig.
From Figure 11, the ratio η of the bell crank actuation to transmission arm has a nonlinear effect on the GII, GVI. The GII increases rapidly to a maximum with the increase of the ratio. The GFI increases almost linearly while the changing trend of GVI is opposite to the ratio. And, the GCPI is insensitive to the ratio change.
Through the above analysis, it can be concluded that the selected geometric parameters have an overall effect on all the performance of the test rig, of which the length of the coupling link and bell crank transmission arm have a significant impact on the coupling of mechanism. As a result, by observing the variation of the established performance indices with an arbitrarily chosen design variable, it can be confirmed that all the four candidate design parameters are critical for the optimal design of the test rig.   therefore defined as where X = [l 1 , l 2 , l 4 , b 11 , η] T is the set of design parameters, λ = [λ 1 , λ 2 , λ 3 , λ 4 ] T is the weighting vector and η = [η h , η f , η v , η cp ] T is the performance indices vector.   The longer the length of the coupling links, the smaller the coupling between the sub-chains, but the stiffness of the material needs to be considered and a compact structure also should be guaranteed. The limits of workspace are selected as W min = [−0.13m, −0.2m, −0.19m] and W max = [0.13m, 0.2m, 0.19m]. The same consideration should be put on the choice of the ratio η along with the size of the transmission arm. As a result, the limits chosen for the design variables are LB = [1.5m, 1.0m, 1.3m, 0.5m, 0.5] T and UB = [1.7m, 1.2m, 1.5m, 0.6m, 0.75] T .The four performance indices mentioned above conflict with each other, as shown in Figures 10-11. Some trade-offs have to be made to select the appropriate structural parameters. In this article, more attention is placed on the decoupling performance index, so the weight vector is selected as λ = [0.1, 0.2, 0.2, 0.5] T .

2) DESIGN RESULTS AND DISCUSSION
The optimization procedure with the ROA algorithm has been carried out in the MATLAB environment, the results of which are presented in the following part.
As can be seen from Figure 12 (a)-(b), the optimal fitness converges quickly and is reached by all particles using about 10 iterations. In Figure 12 (b), the design variables eventually converge to [100mm, 100mm, 100mm, 50mm, 0.5].
The structural parameters before and after optimization by the ROA algorithm are shown in TABLE 8.
In TABLE 8, it can be seen that the length of the coupling links and the bell crank transmission arm all reach the maximum value, while the ratio η reaches the minimum value. To clearly illustrate the change of performance indices before and after optimization, the distribution of the later three indices in the work plane with Z = 0 is shown in Figures 13-15.
In Figures 13-15, it can be observed that the variation of the performance indices is almost symmetrical. The comparison of performance indices before and after optimization is shown in TABLE 9.
After optimization, the speed transmission index κ v increases while the force transmission index κ f shows the opposite trend. The melioration of the speed performance can allow the designer to select the servo valve with the lower-rated flow as well as the oil source with smaller supply flow, which can slash the cost enormously. Although the force transmission performance worsens, the potential of high load capacity of hydraulic systems can be fully utilized. The coupling degree of sub-chains is significantly reduced by optimization, which makes the device more conducive to control.
It can be seen from TABLE 9 that the standard deviations all decrease so the performance becomes more uniform after optimization.

VI. CONCLUSION
A novel root optimization algorithm (ROA) is proposed in this article. Special growing characteristics of roots such as branching, tropism motivates the development of the new algorithm. The branching number is determined dynamically and levy flight enhances the exploration and exploitation capability. The effectiveness of ROA has been proven through comparison with the MFO, GWO, WOA, IWO, GPSO and GA algorithms on classical and CEC 2019 benchmarks. ROA outperforms the mentioned algorithms in several benchmark functions and the statistical results further highlight that ROA is a promising and competitive algorithm. Resorting to the ROA algorithm, a set of optimized dimensional parameters of a multi-axial test rig is obtained successfully. The simulation results show that the test rig achieves an ideal kinematic performance throughout the assigned workspace. The physical prototype of the multi-axial road test rig for further research has been constructed as shown in Figure 3 and actual tests will be carried to validate the performance of the optimized test rig. APPENDIX See Table 10