Proposing a Graphic Simulator for an Upper Limb Exoskeleton Robot

In this study, a graphic simulator that is used to simulate problems related to kinematics and dynamics for an exoskeleton robot arm with 5 degrees of freedom (DoF) was presented. The graphic simulator utilized the advantages of design software SolidWorks, Catia, and the computing and simulation power of SimMechanics Toolbox in Matlab. The core of the proposed graphic simulator is algorithm to solve the kinematics and dynamic problems of a developing upper-limb rehabilitation robot. The study used the proposed optimization-based algorithm to solve the inverse kinematics (IK) problem for the redundant robot model. Endpoint trajectories were imported from measurement data. The joints variable solutions obtained before entering the dynamics problem were smoothed to ensure feasibility in the later calculation process. A process to solve the inverse dynamics problem using physical model by combining the power of two software SolidWorks and SimMechanics was also proposed. This process ensured that the Robot ’ s design could be changed and updated to the kinematics calculation fast and easily. To evaluate this procedure, we also compared these dynamics results with results when applying the Lagrange – Euler formulation. All these calculation and simulation processes have been integrated into the graphic simulator software to show ef ﬁ ciency and user-friendliness.


Introduction
In recent years, with advancements in the fields of mechanics, electricity, control, etc., these are helping to speed up the process of design and production robots in general and rehabilitation robots in particular. Among upper-limb rehabilitation robots, exoskeleton robots are more commonly designed and used than end-effector type robots due to advantages such as rich exercises, precise control of each individual joint, and convenience for measuring and evaluating the effectiveness of the training process [1][2][3]. The main application of upper-limb exoskeleton robots is physical assistance for the rehabilitation process [4][5][6][7][8]. This technology is becoming a very important solution for weak or disabled people in the activity of daily living (ADL). It was developed to enhance the performance and strength of the wear, e.g., Kim et al. [9] and Guidali et al. [10].
Although it has received a lot of attention in recent years in many research, rehabilitation robot design, control, and human-robot interaction aspect are still challenging. Mechanical design and kinematics are the most important issues in developing an ergonomic exoskeleton robot. Meanwhile, dynamic problem is a quite crucial problem in choosing motor and control process. Regarding the inverse kinematics problem, there are some traditional methods such as the geometric method [11,12] and the iterative method. However, when applying these methods to a redundant robot, it often requires complicated and time-consuming calculations. The difficult point is the polynomial as well as the necessary variation of formulas corresponding to different robot structure [13][14][15]. For the upper-limb exoskeleton robot which normally has the same redundancy as the human arm, natural movement of the arm as well as the joints is important. Kim and Rosen [16] suggested the Swivel method which is similar to the geometric method. However, this method first needs to determine the position and orientation of the elbow joint. In addition to the methods mentioned above, in recent years, application of metaheuristic optimization algorithms has become increasingly common [17][18][19]. Recently, Nguyen et al. [21] proposed an algorithm based on the optimal algorithm to solve the inverse kinematics requirement for some human's movements of ADL activities as well as sports activities. In this study, when building graphic simulator software, the author also uses Pro_ISADE algorithm to solve the requirement of inverse kinematics for UExosVN robot.
For the inverse dynamics requirement of robot in general and upper-limb exoskeleton robot in particular, it is commonly solved by Lagrange method [22]. These were mathematical dynamic models. The input information for the method are kinematics parameters, inertial matrix, centripetal and Coriolis torques, and Jacobian matrixes. For complex robotic systems, the methods will need a very large computational volume. Furthermore, when changing robot configuration, the process of updating the parameters is very timeconsuming. Arora and Singh [23] suggested applying physical modeling to model and simulate for a robot. However, the study did not mention whether the material properties when exported from the design model in SolidWorks software were preserved or not. This can affect the accuracy of the results of the dynamics problem. Cong et al. [24] and Vantilt et al. [25] proposed physical modeling parameters getting from the 3D design environment to compare with experimental parameters. The software automatic dynamic analysis of mechanical systems (ADAMS) [26] was used to simulate the dynamics of an upper-limb rehabilitation robot. The results were used to select actuators. In this study, authors proposed a process for calculating dynamics by physical method based on promoting the design power of CAD software and the simulation and calculation power of SimMechanics Toolbox software of Matlab. The strength of the proposed process is to create a connection and unity between kinematics conventions such as zero origin, direction of motion and the process of calculating inverse kinematics, and inverse dynamics for complex models. A computational model is a design model made up of many parts with different materials. Adjustment of design parameters can easily be updated for this calculation to find the required dynamic values.
Relating to the graphic simulator, there are a number of authors who have modeled their robotic system [27,28] to serve for manipulating and testing the control system within the virtual environment. Kereluk and Emami [29] presented a new modular that is autonomously reconfigurable manipulator consisting of 18 DoFs. In this environment, the user can select and change the robot configuration and perform simulation of operation in the virtual environment. López-Franco et al. [31] built a simulation software in Matlab software to simulate the process of applying metaheuristic algorithms to solve the inverse kinematics problem for the robot model. Currently, there is no research to build graphic simulator for rehabilitation robots. The proposed simulators also do not fully handle the problems of kinematics, dynamics, and are less flexible when it comes to mechanical design adjustments.
In this study, the authors proposed to build a graphic simulator that synthesizes kinematics and dynamics problems for a developing upper-limb exoskeleton robot (UExosVN). The robot model consists of 5 active DoFs and 2 passive DoFs. The software combines the strengthen points of 3D design software such as SolidWorks and Catia with the computing power from Matlab, in which, 3D CAD design software provides geometrical and physical parameters such as, mass, center of mass (CoM); moment of inertia, and product of inertia. This information is used to update the physical properties of the model in SimMechanics software to perform dynamic calculation. In the proposed simulators, the Pro_ISADE algorithm was also used to run on the Matlab platform and to calculate the inverse kinematics problem for UExosVN. Following is the structure of remaining paper: The UExosVN model is described in Section 2. Section 3 presents the inverse kinematics and dynamics solutions for the robot. These are the core algorithms inside the proposed graphic simulator software. Section 4 presents the graphic simulator interface. Section 5 presents the results when applying the algorithm and proposed dynamics analysis process. Finally, Section 6 outlines the conclusions.

Forward Kinematics of the
Developing UExosVN 2.1. Structure of UExosVN. The structure of the UExosVN robot was described in [30]. The robot has 7 active and 2 passive DoFs, as shown in Figure 1(a ðaÞ ðbÞ The completed design of the UExosVN model is shown in Figure 1(b).

Forward
Kinematics. The coordinate systems and kinematics parameters of the robots are described as in Figure 2 and Tables 1 and 2, respectively. From that, the D-H parameters of the model can be obtained and summarized in Table 3. The forward kinematics requirement for the manipulator is solved by using Equation (1): where S and C denote the sine and cosine functions. The position and orientation of the end-effector can be determined by Equation (1): The position of end-effector coordinates in the working space is determined by:

ISADE Algorithms with Searching Space Improvement.
Bui et al. [20] and Nguyen et al. [21] proposed the Pro_ ISADE algorithm to handle the inverse kinematics problem for the human arm in daily and sports activities. From this study, it also showed the feasibility of the proposed algorithm in solving the inverse kinematics problem for the exoskeleton robot of UExosVN. A summary of the ISADE algorithm is given in Figure 3. Nguyen et al. [21] proposed a solution to ensure continuity of joints values is to limit the initialization domain of X. The proposed algorithm helps to reduce the calculation burden and ensure the continuity of the joint's values. The target function, in this case, is described in Equation (4) [33].
Applied Bionics and Biomechanics 3 The target function includes two parts. The first one is the difference between the original and the first joints' values. The remaining is squared deviation between the calculated and derived end-effector coordination.
a;b are penalty coefficients. The detail the ISADE and the proposed Pro_ISADE algorithms are shown in Figure 3 and Ref. [21], respectively.
Values q k 1 and q k 0 (i = 1) are the joints' variable values at the original position and 1st point on the trajectory, respectively; x i ; ð y i ;z i Þ and x ei ; ð y ei ;z ei Þ are the end-effector coordinates for the i-point found by the algorithm and the desired end-effector coordinates.

Proposed Process to Solve the Inverse Dynamic for UExosVN.
In this section, the research proposes a process of using physical simulations applied for complex manipulator system [30]. The process has to ensure some following requirements: The proposed process is shown in Figure 4. The process sequence is described as follows: A detail explanation about the requirement and how to proceed for each step is explained in [30].

Evaluate the Proposed Dynamical Model Process.
The results and effectiveness of the proposed algorithm to solve inverse kinematics for the robot were presented in the previous study [21]. In this part, the study will compare the accuracy and advantages of the proposed physical model process with the mathematical method.
where M is the inertial matrix, C is the centrifugal and Coriolis forces matrix, τ is the control torque, and G is the gravitational forces vector; q,q, andq are the vectors of the joint angle, speed, and acceleration joint vectors, respectively. An important point when applying this mathematical method is the method to find the MP for the links. Because the links are composed of parts with different materials, the main moment of inertia values (Ixx, Iyy, Izz) are directed not in the directions of the D-H coordinate axes. The solution to this problem is to use Catia's strengths. In Catia software, it allows users to add a CS (the D-H CS), and then measure the MPs relative to this CS.

Comparison of Physical Modeling Methods and
Mathematical Modeling in Terms of Accuracy. In this study, the author simulated and calculated the inverse dynamics for the UExosVN robot in two complex exercises: drinking water and playing the ball tasks. The reason for choosing these two activities is because they are both complex exercises in daily life: an exercise from activities of daily living and another exercise in sports.
The angle measurement data were performed by the proposed E-HMCS device [21]. The disadvantage of these measurement data is the few numbers of points (about 50 points per trajectory); discrete points are not suitable for performing derivative operations; Measurement data are not compatible with robot arm data.
To overcome the above limitations, the study proposed to perform the dynamic calculation process for the UExosVN robot as shown in Figure 5. This is also Mode 1 in the structure of the graphic interface that will be described in Appendix A. The measured angle data q m i and dq m i i ¼ 1÷5 ð Þ have been smoothed by function spaps x m ; ð y m ;tolÞ in Matlab, where x m ; ð y m Þ is the measured data; tol is tolerable. Tolerances for smoothing q m i and dq m i are suggested as shown in Equation (6), where tol is the tolerance that is set with slider in the graphic simulator that will be presented in later section. dq ms1 i is the velocity after smoothing q m i . The Convert data step is to convert the measuring device data into UExosVN robotic arm data. This kinematics data are fed into the inverse dynamic model to determine the control moment values T R i at the joints.
tol q ¼ −tol=1;000 (1) Case Study 1: Drinking Water Task. Drinking water task was carried out in two phases and depicted in the literature [21]. There are reaching and lifting phases. In the first one, a person moves his/her hand from an initial point to a current point on a table to get a cup. Then, he/she will take the cup and lift it to his/her mouth in the second stage of lifting. When doing exercise all the angle data are collected from the E-HMCS measuring device. Then, the data are smoothed and converted into the robot arm joint variable angles. The values include five joints values which are explained in Table 2. The results about kinematics parameters, error comparison of torque results for different joints by applying the traditional calculation, and the proposed process were depicted before in [30].
(2) Case Study 2: Playing the Ball Task. In this exercise, people moved the racket from a point below of his body diagonally above the body. The kinematics joints result of the UExosVN robot is shown in Figure 6. There is one common thing between drinking water and playing ball tasks is that the joint value 5 does not change. This is explained because the study only focused on the endpoint position while the joint 5 of Elb. Pro./Sup. mainly affects the direction of the endpoint. Similar to the drinking water task, in the playing the ball task kinematics values, as shown in Figure 6, were inputted into the physical modeling method described in Section 3.2 and the mathematical method as described in Section 3.3.1. The results are compared and shown as Figure 7. The results once again demonstrate the equivalence in the application of physical methods and mathematical methods in solving dynamic problems for UExosVN robots.
From the above simulations, the study shows the similarities between the physical simulation methods and the mathematical method when solving dynamic problems for complex robots. The study confirms that dynamics problem can be handled effectively by physical simulation method. The advantages of the method are:

Proposed Graphic Simulator
An operation flowchart of the graphic simulator is depicted in Figure 8. First, we need to choose Mode 1 or Mode 2 of the program depending on whether the available data are information about the joint parameter or endpoint trajectory of the robot, respectively. After that, the kinematics and dynamics of the robot model are processed based on the proposed process as shown in Figures 4, 5, and 9. The physical model for the UExosVN robot is shown as Figures 10 and 11. Another noteworthy point is that, in the design model of the UExosVN robot as described in Section 2.1, the system has 2 passive DoFs which are used to adjust the arm and forearm length to suit each patient. The proposed physical model added these 2 passive DoFs in parameters d 3 and d 5 , as shown in Figure 11. The description about these 5 actives and 2 passive joints as well as those parameters are shown in Section 2.1. The adjustment values are imputed on the graphic interface. Then, the proposed physical model will automatically modify the length of arm and forearm based on the provided values. It is considered as the reconfiguration property. This is also a strong advantage of physical modeling over mathematical modeling.
The operation modes and operation procedure of the proposed graphic simulator software are shown in Appendices A and B, respectively.

Operation
Results of the Software 5.1. Mode 1. After imputing the forearm and arm length parameters, we continue to transfer the joints values and perform tolerance adjustment to smooth. When these parameters are accepted, we press the Calculate button, the kinetic and dynamic results for Mode 1 as described and analyzed in Appendix A are shown.   Figures 12 and 13, respectively. These results were taken to represent one of the five joints of the UExosVN robot arm. These results included the inputted joints value. These were discrete data with about 50 points and quite strong variation. The smoothing process created a continuously joint variable with a fixed sampling period. It facilitates the computation of subsequent velocity and acceleration parameters. The kinematics and control moment results for each link were also shown. The graphic simulator also contained a graph comparing the endpoint trajectories before after smoothing, as shown in Figures 12(b) and 13(b).

Case Study 2:
Playing the Ball Task. Performing completely similar to the case of drinking water task, we get the results of kinematics, dynamics, and trajectory of the endpoint of the task as shown in Figure 14.

Mode 2.
The calculation procedure in this case has been described in detail in Appendix A and Figure 9. In this process, the software adds the solution of the inverse kinematics problem by the Pro_ISADE algorithm. The results obtained were the joints variable's values and then they are also processed with the same basic steps as Mode 1. Inverse dynamic  Torque   T5   T4   T3   T2   T1   T1   T2   T3   q1   dq1   ddq1   q2   dq2   ddq2   q3   dq3   ddq3   q4   dq4   ddq4   q5   dq5   ddq5   T4   T5  T5   T3   T1   T2 FIGURE 10: Physical model for the UExosVN [30]. The results in this exercise include the values of the joint variables from the IK problem, the kinematics values, and the dynamic values of each joint. Each of these results for reaching stage and lifting stage is depicted in Figures 15 and 16, respectively.

Joint parameters
In the results of the comparison about endpoint trajectories, as shown in Figures 15(b) and 16(b), the software presented three trajectories corresponding to the inputted raw data trajectories, the smoothed trajectories, and the trajectories after smoothing joints values.

Case Study 2:
Playing the Ball Task. Performing completely similar to the case of drinking water task, we get the results of kinematics, dynamics, and trajectory of the endpoint of the exercise as shown in Figure 17.

Discussion and Conclusion
The research has been conducted to build a graphic simulator for calculating kinematics and dynamics for an exoskeleton rehabilitation robot arm consisting of 5 active DoFs and 2 passive DoFs. The software has included components of the calculation process such as forward kinematics, inverse kinematics, and inverse dynamic. To solve the inverse kinematics problem, the study used the Pro_ISADE algorithm. This algorithm has been proven by the author to be effective for robot UExosVN. In addition, the study also proposed a procedure to calculate the inverse dynamics problem by physical modeling method applied to complete and complex robot design models containing many details made from many different materials. The physical model allows for quick adjustment and flexible configuration changes. Simulation results for two exercises including daily activities and sports activities demonstrated the effectiveness of the proposed procedure. In addition, when applying these processes to the software, it also provides a user-friendly, intuitive interface that is easy to adjust and get the full results of the kinematics and dynamics of the robot arm such as joints values, joints velocity, and joints acceleration; torque controls the joints; robot arm endpoint trajectory.
The main weak point of the current software is the path planning for the endpoint of the robot. As aforementioned, the   Applied Bionics and Biomechanics 13 data of end-effector trajectory were recorded and imputed to the software by using the human motion capture system. To improve the software, the next research will focus on solving the problem of generating the trajectory for the robot's endpoint. Several solutions can be used like learning by demonstration, reenforcement learning, etc. After solving the problem of setting the endpoint trajectory, this module will be added to the current software to improve the graphic simulator software.

A. Operation Mode of the Software
The graphic simulation interface for UExosVN's calculations is depicted in Figure 18. The software has two operation Modes: (i) Mode 1: Input values are joints' values. The processing sequence has been described in Section 3.3.2. (ii) Mode 2: Input is the endpoint trajectory. The processing sequence is described, as shown in Figure 9. Unlike Mode 1, in this case, raw data about endpoint coordinates x m e ; ð y m e ;z m e Þ are smoothed first to achieve data x m es ; ð y m es ;z m es Þ. These data are further fed through the inverse kinematics problem and processed by the Pro_ISADE algorithm which is described in Section 3.1 to obtain raw joints' variables parameters q R i i ¼ 1÷5 ð Þ . These data are further fed into the smoothing process as well as solving the inverse dynamics problem by physical modeling method to achieve the necessary dynamic parameter T R i .

B. Operation Procedure of the Software
First, we need to determine the length of the forearm and arm segments and input them into the graphic simulator as in Figure 18. Next, depending on the type of information we have, is the value of the joints' variables or endpoint trajectory of robot UExosVN, we choose the operating mode as Mode 1 or Mode 2, respectively.
(i) In the case of selecting Mode 1: The procedure is similar to Figure 5. One note is that the smoothing phase of the trajectory, velocity, is done by moving the slider of Adjust Tolerance and then pressing the Check button to update whether the smoothed level on the corresponding graphs meets the requirements. (ii) In case of selecting Mode 2: The procedure is done as shown in Figure 9. The difference from Mode 1 is that in this case, the software needs to smooth the trajectory and calculate the inverse kinematics. In addition, this process does not require a conversion process.

Data Availability
Data supporting this research article are available from the corresponding author on reasonable request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

Funding
The research reported in this paper was supported by the Ministry of Science and Technology of Vietnam, under award number (ĐTĐLCN. 28/20). This work was also supported by the Centennial SIT Action for the 100th anniversary of Shibaura Institute of Technology entering the top 10 at the Asian Institute of Technology.