Abstract

This paper presents a novel redundantly actuated 2RPU-2SPR parallel manipulator that can be employed to form a five-axis hybrid kinematic machine tool for large heterogeneous complex structural component machining in aerospace field. On the contrary to series manipulators, the parallel manipulator has the potential merits of high stiffness, high speed, excellent dynamic performance, and complicated surface processing capability. First, by resorting to the screw theory, the degree of freedom of the proposed parallel manipulator is briefly addressed with general configuration and verified by Grübler-Kutzbach (G-K) criteria as well. Next, the inverse kinematics solution for such manipulator is deduced in detail; however, the forward kinematics is mathematically intractable. To deal with such problem, the forward kinematics is solved by means of three back propagation (BP) neural network optimization strategies. The remarkable simulation results of the parallel manipulator demonstrate that the BP neural network with position compensation is an appropriate method for solving the forward kinematics because of its various advantages, such as high efficiency and high convergence ratios. Simultaneously, workspaces, including joint space and workspace of the proposed parallel manipulator, are graphically depicted based on the previous research, which illustrate that the proposed manipulator is a good candidate for engineering practical application.

1. Introduction

Various types of mechanisms have been applied in many robotic fields. Nowadays, some parallel manipulators have attracted exhaustive attention from both academia and industry. Compared to series manipulators, parallel manipulators are claimed to possess inherent advantages such as high stiffness, high loading capability, high precision, low error accumulation, fast response speed, and high orientation capability. These characteristics allow the parallel manipulator to be widely used in various alternative fields, such as high-speed machine tools, high-speed pick and place applications, surgical manipulators, and force/torque sensors [13]. However, parallel manipulators suffer inherently from an unfavorable workspace. Therefore, it is of crucial significance to make full use of the potential parallel manipulators to satisfy the current work requirements [4, 5].

Generally, some performance indices, such as dexterity, stiffness, accuracy, and workspace, will be involved in the design process of the parallel manipulator. Kinematics solutions, establishing the relationship between the joint space and operation workspace including position, velocity, and acceleration, are an extremely key issue. The inverse kinematics solution is to find the joint angles or actuator lengths when the position and orientation of the moving platform are given [6], whereas the forward kinematics solution is to find the terminal position and orientation in consideration of a given set of joint angles or actuator lengths. Due to the inverse kinematics, a mathematical solution can be obtained directly through the closed-loop vector of the parallel manipulator, but the forward kinematics solution is very complex and nonlinear, which presents a considerably difficult problem. Therefore, computing the forward kinematics rapidly is critically crucial in real-time control of the redundantly actuated 2RPU-2SPR parallel manipulator. In view of the literature, there is an abundance of research on the forward kinematics. Li and Xu [7] adopted the Newton iterative method to solve the forward position kinematics of the 3-PRS parallel manipulator. Fu et al. [8] studied the polynomial solutions of the forward kinematics problem, converted it into a first-degree polynomial equation and an eight-degree polynomial equation, and verified the effectiveness and correctness of the proposed algorithm. Joshi and Tsai [9] proposed a methodology for solving the forward kinematics of a class of 3-DOF, 4-legged parallel manipulators in closed form and demonstrated the algorithm with the Tricept manipulator that has at most twenty-four real solutions. Wu et al. [10] presented a simple numerical method for forward kinematics of a general 6-DOF parallel, which can generate a unique actual solution. The method possesses more advantage than the conventional numerical iteration method does because it does not rely on the initial value. Merlet [11] presented an efficient alternative algorithm, interval analysis, which can be numerically solved for the forward kinematics, that is, to determine all possible poses for given joint variables. Tsai et al. [12] employed Bezou’s elimination method and optimization techniques to solve the forward kinematics problem of the 3-PRS parallel mechanism, which makes it possible for real-time control applications. Karimi and Nategh [13] utilized a statistical approach including Bates and Watts to investigate the nonlinearity of the forward kinematics, and the results demonstrate that the length of the region has a significant impact on the nonlinearity of the parallel manipulator. Li et al. [14] proposed a novel overconstrained three-degree-of-freedom spatial parallel manipulator and derived the analytical formulation of the forward kinematics. Wang et al. [15] proposed a novel forward kinematics algorithm, which is verity in the 3-RPS parallel mechanism, and the accuracy is, less than 10-6, sufficiently high. However, it needs to select a reasonable initial value.

Recently, neural networks are employed to solve the forward kinematics of a parallel manipulator and draw particular interests for numerous researchers, due to their considerable ability to approximate nonlinear mapping functions. Parikh and Lam [16] proposed several neural network strategies for the Gough-Stewart parallel manipulator, and results show that the hybrid strategy, combined with a standard Newton-Raphson numerical method, is the best methodology. Amir et al. [17] proposed a wavelet neural network, which can be employed as a powerful method for nonlinear mapping equations, to approximate the desired trajectory of a hybrid robot manipulator. Sadjadian et al. [18] used neural network identification with different structures to solve the forward kinematics considering highly coupled and nonlinear characteristics of the mathematical model which are almost impossible to solve analytically. Cheng et al. [19] generate the workspace of the 3-PRS parallel robot by adopting the trained BP neural network approach. However, it is rare to use the neural network for solving the forward kinematics of the redundantly actuated parallel manipulator.

The goal of this paper is to solve the forward kinematics, which promotes further practical applications of the promising redundantly actuated 2RPU-2SPR parallel manipulator. The remainder of the paper is organized as follows. The research background of the forward kinematics solution is presented firstly. In Section 2, a novel redundantly actuated 2RPU-2SPR parallel manipulator is proposed, which can be selected as the main body of the hybrid machine tool. In Section 3, the degree of freedom of the proposed parallel manipulator was briefly analyzed in a general configuration via the screw theory, and the modified Grübler-Kutzbach (G-K) equation is utilized to verify the correctness for the degree of freedom. Then, the inverse kinematics position is formulated. The BP neural network optimization strategies proposed here are trained to solve the forward kinematics based on the inverse solution in Section 4. In Section 5, the joint space and workspace of the parallel manipulator are performed, respectively, and the workspaces of the parallel manipulator are intuitively depicted by using MATLAB programming. Finally, some concluding remarks are provided in the last section.

2. Manipulator Description

A five-axis hybrid kinematic machine tool is design to mill a large complex structure component surface in the aerospace field. As shown in Figure 1, it is composed of a novel three-degree-of-freedom parallel manipulator and two adjustable annular guides. The three-degree-of-freedom parallel manipulator is a novel redundantly actuated parallel manipulator, which consists of a moving platform connecting a fixed base by four identical driving legs (or limbs) with a symmetrical structure. Leg 1 and leg 3 connect by means of a fixed base to the moving platform by a revolute joint (R), a prismatic joint (P), and a universal joint (U) in sequence forming an RPU configuration kinematic chain, while leg 2 and leg 4 connect by means of a fixed base to the moving platform by a spherical (S), prismatic (P), and revolute joint (R) in sequence forming an SPR configuration kinematic chain, where the prismatic joint is activated independently by a linear serve motor. The kinematic model of the parallel manipulator is shown in Figure 2. The moving platform and the base are a square layout, and its circumcircle radii are and , respectively. The kinematic joint configuration arrangement satisfied certain geometric constraints, and the moving platform will own three degrees of freedom, namely, one translation and two rotations [20, 21].

As shown in Figure 2, the Cartesian fixed coordinate reference frame is attached to the centroid of the fixed base, whose coordinate origin is a square center of the circle, in which the and axes are the same as the vector and , and the axis is determined by the right hand rule. Analogously, the local coordinate reference frame is considered at the centroid of the moving platform. The coordinate origin is the geometric center of the moving platform, the and axes are the same as the vector and , and the axis is determined by the right hand rule. The position vector of point with regard to the Cartesian fixed coordinate reference frame can be indicated as ().

3. Kinematics Analyses

In this section, the kinematics analysis, mobility, inverse kinematics, and forward kinematics of the redundantly actuated 2RPU-2SPR parallel manipulator are mentioned.

3.1. Mobility Analysis

The mobility determination is the first and foremost issue in designing the motion characteristic of a parallel manipulator. In order to determine the motion pattern of the 2RPU-2SPR parallel manipulator, mobility analysis including the number of degrees of freedom and motion classification is indispensable. Furthermore, it is necessary to analyze the constraint wrench screw provided by each leg to the moving platform and comprehensively analyze the constraint type of the moving platform, i.e., whether they are a constraint force or constraint moment. Therefore, the detailed analysis process is as follows [22].

In Figure 2, through the analysis of the screw theory, the twist screw of the first RPU leg can be expressed as follows in the fixed coordinate system where and represent the direction cosine of the drive joint, and and represent the direction cosine of the third rotation axis of the universal joint of the 1-st leg.

Then, considering the reciprocal product, the wrench system of equation (1) can be easily obtained as where presents a constraint force that passes through point and is parallel to , and represents a constraint moment perpendicular to and .

Analogously, the twist screw of the second SPR leg can be expressed mathematically as

The wrench screw of equation (3) can be described as where represents a constraint force that passes through point and is parallel to .

The twist screw of the third RPU leg can be expressed numerically as

The wrench screw of equation (5) can be obtained as where represents a constraint force that passes through point and is parallel to , and represents a constraint couple perpendicular to and .

In a similar approach, the twist screw of the fourth SPR leg can be expressed mathematically as

The wrench screw of equation (7) can be obtained as where represents a constraint force that passes through point and is parallel to .

So far, the constraint wrench screw of the redundantly actuated parallel manipulator can be expressed as

It is noteworthy that the direction vector coincides with the revolution joint axis of the moving platform, whose information plays an important role in this additional potential relation, which can be described as

By combining equation (10) and rearranging equation (9), here the twist screw of the moving platform can be formulated as follows:

In equation (11), represents one rotational degree of freedom around the x axis, represents a translational degree of freedom that is perpendicular to the moving platform, and represents one rotational degree of freedom of the moving platform around the axis.

Generally, the number of degrees of freedom of the parallel manipulator can be calculated by utilizing the modified Grübler-Kutzbach (G-K) criterion, that is, where represents the degree of freedom of the mechanism, represents the number of the components, represents the number of the kinematic joints, represents the order of the mechanism, represents the degree of freedom of the -th kinematic joint, represents the redundant constraints of the mechanism, and represents the local degree of freedom.

There is neither a constraint couple in the same direction nor a constraint force in collinearity among all the wrench screws of the parallel manipulator; therefore, there is no common constraint, that is, . Because there are only three linearly independent variables in the six constraint screws of the parallel manipulator, therefore, the parallel manipulator has three redundant constraints, that is, . As the parallel manipulator has no local degree of freedom, so . We can see from the schematic of the mechanism that the number of the component is 10, the number of the kinematic joint is 12, and the relative freedom of all the kinematic joints in the mechanism is 18.

Thus, the mobility number of the 2RPU-2SPR parallel manipulator can be determined easily by utilizing equation (12):

In summary, the redundantly actuated 2RPU-2SPR parallel manipulator has three independent degrees of freedom, that is, two rotational degrees of freedom around the axis and axis and one translational degree of freedom along the axis.

Compared with the existing 4-DOF parallel manipulator possessing four actuators [23], this novel redundantly actuated parallel manipulator has some key merits such as four actuator joint variables with redundant actuation and three redundant constraints, which make us draw more attention to its advantages, i.e., increase the desirable workspace, improve the dexterity, and eliminate/decrease the singularity. Here, each universal joint is composed of two mutually perpendicular rotational joints with “T” configuration type not “+” configuration type, and the “T” configuration type universal joints designed in the mechanism can present a relative greater rotation angle which can avoid singularity, increase the workspace, and improve the kinematics and dynamics performance [24].

3.2. Inverse Kinematics

Figure 2 shows the structural diagram of the parallel manipulator, where the closed vector constraint equation can be written with regard to the fixed coordinate system

In equation (14), is the position vector in the fixed coordinate system, represents the linear displacement of the -th actuator, is the position vector in the relative coordinate system, and is the representation of vector with respect to the fixed coordinate system, namely, , where is the rotation matrix of the parallel manipulator. From the structure of the parallel manipulator, the position vector can be written as

As illustrated in Figure 2, the homogeneous coordinate transformation matrix of the redundantly actuated parallel manipulator can be obtained through three times of complex transformation, and the moving process is as follows. First, angle is rotated around the axis. Second, distance is translated along the new axis. Finally, angle is rotated around the axis of the moving platform. Its pose can be represented by position vector and rotation matrix , and the transformation matrix can be expressed as

It can be seen from equation (16) that the position vector of the moving platform is a function of independent variables , , and , namely,

So, the parasitic motion in the movement of the parallel manipulator is related to the workspace parameters, that is,

Equation (18) represents the coupling relationships between the position and orientation variables. The three workspace variables , , and are independent of one another, so the inverse kinematics solution can be calculated with the following analytical equation: where represents the Euclidean norm of a vector, so the unit vector in the direction of activate leg can be expressed as

Taking the positive square root of the resulting equation (19) and rearranging it lead to

If the structure parameters of the parallel manipulator and pose of the moving platform are given, the displacement of the activated joint can be calculated by employing the inverse kinematics equation (21). The purpose of the inverse kinematics solution is to determine the actuator lengths when the position and orientation (pose) of the moving platform are given. It is such length and pose data that are adopted as the training samples of BP neural networks in the next section.

4. Forward Kinematics of the BP Neural Network

The forward kinematics solution is aimed at determining the position and orientation of the moving platform with respect to the fixed reference frame once the linear displacement of the joint variables is given. In consideration of a rotation axis of the proposed parallel manipulator located at the axis around the fixed coordinate system, there is parasitic motion in the movement process, which makes the forward kinematics solution of the parallel manipulator especially complex. Therefore, the forward kinematics solution needs to be calculated by applying an artificial intelligence method named the BP neural network in this section [25, 26].

4.1. BP Neural Network

The BP neural network is a dynamic system of the topological structure of the digraph, which employs differentiable function as activation function to propagate forward. At the same time, the weight factor (, and presents the node from the input layer to the hidden layer as well as from the hidden layer to the output layer) of network parameters is adjusted continuously through error back propagation to minimize the mean square errors of the network and realize the nonlinear mapping from the actuator joint variable space to the operation variable space without revealing the mathematical equation beforehand. Its topological structure, as shown in Figure 3, contains the input layer, hidden layer, and output layer. In this case, the input layer consists of four nodes, i.e., four linear input actuator displacements , and the output layer consists of three nodes, i.e., one position and two orientations parameter [27, 28].

The number of the nodes () in the hidden layer is generally determined by empirical formulas [29], that is,

Among them, is the number of nodes in the input layer, is the number in the output layer, and is the constant whose range is from one to ten. The number of nodes in the hidden layer is selected as six in this paper.

The differentiable and linear functions are usually employed as alternative activation function of the BP neural network. In this paper, the S-type logarithmic function, Tan-Sigmoid, is selected as the activation function of output hidden layer nodes (neurons) in the prediction model, as follows:

Back propagation requires that the activation function used by the artificial nodes be differentiable, so the derivative of the Tan-Sigmoid function can be expressed as

Since the magnitude of different data is not the same, which may result in slower convergence ratios and long training time of the neural network, furthermore, the range of activation function in the output layer of the neural network is limited to one interval, so it is necessary to normalize the target data of the training network to the range [-1, 1]; the normalization criterion can be written as where is the minimum value of , is the maximum value of , is the input vector, and is the output vector of normalization.

The learning and training samples are set up by utilizing the inverse solution equation (21) of the parallel manipulator. In this paper, we choose 5512 data values representing the learning samples; the network training function is traingdx which updated weight and bias values according to the gradient descent momentum and an adaptive learning rate. The network performance function is the mean square error, which is the error between the network calculated outputs and the target outputs. The number of iteration epochs is set as 4000, the expected error goal is , the learning rate is 0.01, and the training network is trained by the Levenberg-Marquardt (LM) optimization algorithm to improve the training speed and reduce the convergence error.

After being successfully trained, once all the kinematic parameters are known, the forward kinematics can be directly computed, which is the main goal of this paper. For the sake of convenience, twelve groups of positions and postures are taken as the test data, their inverse solutions are brought into the trained network, taking the first group set of data as an example, and the iterative process is depicted in Figure 4. After 722 epochs, the optimum parameter values () can be obtained within the allowable error range, which implement the nonlinear mapping from the actuator joint variable space to the operation variable space [30, 31].

The forward kinematics solution of the parallel manipulator can be solved by minimizing the network performance function. The mean square error (MSE) of the BP neural network can be derived according to the error back propagation equation where is the predicted value of the BP neural network, and is the desired target value; here represents parameters , , and .

All of the selected parameters (, , and ) and the errors calculated by the BP neural network are shown in Table 1.

From Table 1, we can see that the range of errors is (, ); the magnitude of most errors is almost in the order of 10-3. However, the error of the third group and the sixth group is smaller than others, their magnitude of error is 10-4, and their orientation angles are the same, which indicates that the mean square error is very small when the orientation angle is in the same value. We also find that in the eleventh group, the error is the biggest of all the results, and their orientation angles are in different directions, so the error is particularly large. By comparing the fifth group and the ninth group, it is illustrated that the mean square error also becomes larger with the increase in the orientation angle. We hope that a smaller error is much better; therefore, much more efforts should be devoted to decrease the mean square error.

4.2. Improved BP Neural Network

Since the position and orientation values of the parallel manipulator are not in an order of magnitude, therefore, the improved BP neural network is used to train position and orientation ()separately, and the LM optimization algorithm is also used in the same manner to solve the forward kinematics problem. The data of the above twelve groups are predicted with the improved BP neural network. The results of the parallel manipulator for calculation are shown in Table 2.

The results of the improved BP neural network are more accurate than those of the traditional BP neural network. The order of magnitude has basically reached 10-4, and there has been a great improvement compared with the traditional BP neural network. However, the accuracy is still difficult to satisfy the high accuracy requirements of the parallel kinematic machine tool. Therefore, error compensation is needed to further improve the calculation accuracy.

4.3. BP Neural Network with Position Compensation

Even though both BP neural network and improved BP neural network approaches yield good results for training and testing for given desired data points, in order to obtain a more satisfactory accuracy, a BP neural network algorithm with position error compensation is adopted; the calculation process in Figure 5 represents the sequence followed for minimizing the mean square error to obtain the optimum forward kinematics.

Step 1. Select the actuator displacement as the initial input value of the loop program.

Step 2. Calculate the pose of the moving platform with the improved BP neural network.

Step 3. Solve the actuator displacement with the inverse kinematics solution.

Step 4. If the delta absolute value less than the permission accuracy , then export the pose parameter ; otherwise, let as the new input of next cycle repeat step 2 until we obtain the satisfactory values.

Analogously, the forward kinematics solution can be derived by employing the BP neural network with position compensation, and the approximation permission accuracy ; after calculation, the forward kinematics solution results and mean square error can be obtained, as shown in Table 3.

From Table 3, we can see that the accuracy of the forward kinematics solution is greatly improved by using the BP neural network with the error position compensation method. The magnitude of the error has already reached 10-6, which indicates that this optimization strategy has higher accuracy than the previous two methods described above. The calculation results are shown in Table 3. It can be seen from the table that the maximum mean square error is , while the minimum one is . The experimental accuracy has been significantly improved, which also illustrates the effectiveness of the excellent algorithm for calculating the forward kinematics.

Compared with three neural network strategies, as shown in Figure 6 (the 1st represents the general (traditional) BP neural network, the 2nd represent the improved BP neural network, and the 3rd represents the BP neural network with position compensation), we found that the BP neural network with position compensation has the best effect for the forward kinematics solution of the parallel manipulator. Comparing the 1st method with the 2nd method, the mean square error of the improved BP neural network method has much more improvement than that of the 1st method.

However, as depicted in Figure 7, the traditional BP neural network has a much lesser mean square error in the direction and a relatively bigger orientation error including angles and , while the improved BP neural network method is just opposite to that of the 1st method. It is also shown that the 1st and 2nd are not good enough, and either method has shortcomings and inadequacies. Comparative analysis illustrated that the BP neural network with position compensation is fairly efficient and possesses high accuracy. It is also proved that the proposed algorithm can be efficiently applied to forward kinematics of the 2RPU-2SPR redundantly actuated parallel manipulator, which is critically important for the further research including dynamic analysis, control simulation, and experiment prototype real-time control.

5. Workspace Analysis

The determination of joint and workspace is one of the most important issues of the redundant actuation parallel manipulator, which directly reflects the comprehensive performance of the parallel manipulator; its analysis results provide a theoretical basis for the design and application of the 2RPU-2SPR parallel manipulator [3234]. A real experiential prototype (as depicted in Figure 8) is fabricated, whose workspace is generated by adopting the main loop algorithm based on the forward kinematics combining the structural parameters and the limitation conditions.

For convenience, we summarize the main loop algorithm procedure as follows.

for to by step
  for to by step
   for to by step
    for to by step
Forward_kinematics = BP_compensation (, , , );
  If iteration number ‖ error accuracy<eps
break;
   end
record (, , , , ,,);
    end
   end
  end
end

In the actual movement process of the parallel manipulator, the following limitations should be considered: (1)Limitation of the actuated joints is defined as where and represent the minimum and maximum actuator displacement, respectively.(2)Joint angle constraintsUnder consideration, the rotation angle of the revolute joint and the rotation angle of the spherical joint should satisfy the following inequality: where , , , and represent the minimum and maximum constraint angles of the revolute joint and the spherical joint, and represents the normal vector of the fixed base.(3)Singularity configuration constraintAt a singularity configuration, the parallel manipulator may lose one or more degrees of freedom, which will seriously affect the kinematic performance of the parallel manipulator, and then those configurations should be avoided. So the joints are kept away from the singularity configuration, and the rank of the constraint screw system is limited to be three, that is,

As has been described beforehand, the architectural parameters and the maximum and minimum limit values of constraint conditions are shown in Table 4 as below, so the joint space and workspace can be depicted by MATLAB programming.

With consideration of four actuators, the joint workspace can be combined by (), (), (), (), and so on; for the sake of convenience, we just provide only one joint workspace () to illustrate the forward kinematics solution, as depicted in Figures 9 and 10, in which there is no empty and has a smooth surface and envelope boundary. Figure 9 indicates the joint space that actuator joints can reach the points set. Figure 10 is the projection in the dimensional plane.

Simultaneously, the boundary of the workspace can be formulated by adopting the BP neural network with position compensation when given the feasible actuator joint variable displacement, which denotes the available utmost range of the end-effecter of the parallel manipulator. Figure 11 shows the operation workspace of the redundantly actuated parallel manipulator, and Figure 12 shows the projection of the two-dimensional view. From the two figures, we can see that the distribution is uniform, the change is stable, and there is no sudden change, which illustrate that the proposed parallel manipulator has a greater workspace and comparatively better orientation operation capability. Therefore, the parallel manipulator can reach to a certain position with one or several different orientations with good dexterity, which is just required to further improve the flexibility in completing the complicated surface machining. What is more, due to the symmetry of the parallel manipulator, the distribution of the operation workspace about the coordinate axis also presents certain symmetry.

6. Conclusions

In this paper, a novel 1T2R redundantly actuated 2RPU-2SPR parallel manipulator has been presented, which can be applied to high-speed machining in a serial and annular guide rail for milling some large heterogeneous complex freedom surfaces in the aerospace field. According to the investigation of this paper, the following contributions can be drawn: (1)The mobility of the proposed redundantly actuated 2RPU-2SPR parallel manipulator has been analyzed in detail based on the screw theory and Grübler-Kutzbach (G-K) criterion, and the correctness of the analysis method is verified(2)The inverse kinematic position of the redundantly actuated parallel manipulator has been sequentially analyzed, which lay the basis for the BP neural network training samples. The forward kinematics solution problem of the redundantly actuated parallel manipulator can be effectively solved by the BP neural network with position compensation with high accuracy, and the calculation results are in good accordance with the high precision requirement of the milling process, which lay a foundation for further dynamic modeling and real-time control of the experimental prototype(3)The joint space and workspace are depicted and constructed. The analytic results illustrate that the proposed manipulator is a good candidate for engineering application. The proposed methodology to realize two rotations and one translation can also extended to other parallel manipulators including redundantly actuated and nonredundantly actuated parallel manipulators with a complicated forward kinematics problem as well. In the future works, we will focus on the dynamic characteristic analysis, error identification and compensation, sensitivity analysis, and real-time control

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Acknowledgments

The authors would like to acknowledge the financial support of the Fundamental Research Funds for the Central Universities under Grant Nos. 2018JBZ007, 2018YJS136, and 2017YJS158; China Scholarship Council (CSC) under Grant No. 201807090079; the Natural Sciences and Engineering Research Council of Canada (NSERC); and the York Research Chairs (YRC) program. Meanwhile author Haiqiang Zhang is grateful to Advanced Robotics and Mechatronics Laboratory and the science librarian John Dupuis in York University.