Musculoskeletal modeling and humanoid control of robots based on human gait data

Human gait data collection

At present, the motion capture system is divided into five categories according to the principle: mechanical motion capture system (Wu et al., 2005), acoustic motion capture system, electromagnetic motion capture system (Guo et al., 2011), inertial motion capture system (Kim & Nussbaum, 2013) and optical motion capture system (Kurihara et al., 2002; Guerra Filho, 2005; Kirk, O’Brien & Forsyth, 2005). Among them, the optical motion capture system is divided into two categories: motion capture system based on computer vision (optical non-calibration) and optical motion capture system based on mark point (optical calibration).

This paper selected Nokov 3D infrared passive optical motion capture system with high accuracy and good effect after comparing several existing motion capture systems and combining with the needs for current research topic. In the scene set up by this system, the infrared camera is used to fully cover the experimental scene, infrared light is emitted by the infrared camera array in the process of data collection, and the position information of the reflective Maker points are captured in the experimental scene. In the process of collecting gait data, the experimental subject walks in the experimental scene with affixed Maker points. In order to meet the needs of the research, the experimental platform is equipped with a three-dimensional force measuring platform, which can synchronously collect the three-dimensional ground reaction force during the movement of the experimental object.

Before collecting experimental data, the deployment of the experimental platform is also critical. The deployment of the camera position has a fatal impact on the experimental data (Kurihara et al., 2002). In the experiment, the influences of different camera arrangement methods on experimental data were tested. It was found that the best data is obtained by using the approximate circular camera arrangement. This arrangement allows each camera to maximize its utilization value. In the experiment, the cameras scene is arranged around the force measuring platform in an oval shape. The calibration origin is positioned as far as possible in the center of each camera’s field of view by adjusting the orientation of the cameras. The experimental collection scene is shown in Fig. 2.

Gait data processing

After the data collection is completed, it is preprocessed to ensure the completeness and accuracy of each frame of data. For missing data points, we had appropriate discarding or interpolation methods for processing. For severely missing data, the entire group will be deleted without applying.

The force was collected by three-dimensional measuring platform that is the force between foot and ground when the experimenter walks. During the process of gait collection, the force in the vertical direction is the most important force, it reflects the phenomenon of overweight and weightlessness during the gait cycle. The three-dimensional force as shown in Fig. 3. It can be found that the force between left foot and right foot with the ground is basically symmetrical on same axis. The data of the force platform is zero before the foot contacts it. Next, a small fluctuation is formed in the coronal and sagittal axes at first, then increases to a maximum and gradually decreases to zero. On the vertical axis, it increases rapidly, forming a bimodal curve similar to M, and then rapidly returns to zero too.

The acceleration of the left and right joint has certain symmetry and periodicity, as shown in Fig. 4, it can be seen from the figure that the acceleration from the left leg joint and the right leg joint can basically coincide with each other in the case of 0.5 s of translation. During the walking process of the subject, the acceleration of hip joint points were maintained from 1 m/s² to 5 m/s². The acceleration of knee joint points were maintained from 0 m/s² to 13 m/s². The acceleration of ankle joint points were maintained from 0 m/s² to 25 m/s². It can be easily observed that the acceleration at the ankle joint points is greater than the acceleration at the knee joint points, and the acceleration at the knee joint points is greater than acceleration at the hip joint points. During walking, consistent with experience, the further away from the torso, the acceleration of the nodes greater. Here, a small idea is proposed: based on the acceleration of the joints, a body movement comfort function is designed to evaluate the patients’ comfort in the process of lower limb exoskeleton rehabilitation training. This will be a research direction in the next stage of this subject.

Musculoskeletal model scaling

In this paper, the open-source musculoskeletal model was scaled to obtain the exclusive model equivalent to the experimental object. In order to ensure the accuracy of model scaling, the model was scaled several times. Finally, the accuracy of the model with a scaling error of less than one-thousandth is achieved.

Model scaling allows the open sources model to match our experimental subjects as closely as possible. In the scaling process, static data collected are mainly used in this paper (the experimental data collected while the experimental object is standing still). Before scaling the model, Maker points were added at the appropriate position on the model in accordance with the experimental object. As shown in Fig. 5. Meanwhile, these Maker points were connected to the model bones. In order to ensure the scaled model more accurate, the collected action data are processed in this paper. Through calculation, the left and right width of the pelvis and the length of the left thigh, left calf, right thigh and right calf were obtained(as shown in Table 2), where the mass and length are calculated through Zatsiorsky (Vaughan, Andrews & Hay, 1982) and Harless study (Drillis, Contini & Bluestein, 1964). Finally, the length of these body segments in the model was built.

All the above body segment lengths were obtained from the processing of experimental data, and the measured body segment lengths were averaged. The comparison shows that the data is relatively accurate. In this scaling process, we preserve mass distribution during scale, and change the scale weight of the makers to get a more accurate model.

Inverse Kinematics (IK)

Forward kinematics calculate the final position of the model by giving the initial position, velocity and acceleration of the model. The IK are opposite to the forward dynamics. IK figures out the motion process of the model based on the given position, including the change of physical quantity such as velocity and acceleration in this process.

The IK uses the motion capture data collected in the experiment (walk.trc). And the internal algorithm of the software was used for biological simulation, and the inverse solution was used to calculate the joint angle, pelvis tilt, et al. Among the lower limb movements of the human body, the hip joint movement is most complicated. The hip joint has three degrees of freedom. Therefore, the leg will perform three axial movements on the hip joint, including flexion and extension of the hip joint on the sagittal plane, adduction and abduction on the coronal plane, and internal rotation and external rotation direction of the thigh. During the gait cycle, the hip joint angles change as shown in Fig. 6. The hip joint of flexion is the movement in the sagittal plane, and its range of change is stable at 0 ∼ 0.8 rad. This movement change is the leading and effective movement of the hip joint during the gait, and the movement in this direction will drive the body to move forward. The hip joint of adduction changes smoothly in the gait, with only a slight fluctuation. The hip joint of rotation includes internal rotation and external rotation of the thigh during the gait. It can be seen that the data in this part has strong characteristics, and the range of variation is stable at −0.6 ∼ 0 rad. The hip joint of adduction and the hip joint of rotation have more personal characteristics related to personal habits and leg health conditions. It is also a key factor in judging whether the gait is abnormal. Through IK, the collected motion capture data can be matched with the calibration data of the experimental object, and the motion simulation process of the model can be optimized. The IK tool calculates universal coordinate values for each time step (frame) of the movement. Then, the model is positioned in the “best match” pose of the experimental marker time step. In other words, mark points in the collected motion process are matched with the motion capture data, and the weighted square error of the mark points and motion capture data is minimized.

The law of weighted least squares problem during IK solved by the function: (1)${\displaystyle \underset{q}{min}\left[\sum _{i\in markers}{w}_i\parallel {x_i}^{exp}-x_i\left(q\right){\parallel }^2+\sum _{j\in unprescribed-coords}{\omega }_j{\left({q_j}^{exp}-q_j\right)}^2\right]}$

Where, q is the vector of the generalized coordinates solved, X_i^exp is the experimental location of mark i, X_i(q) is the position of the corresponding mark points in the model (depending on the coordinate value), q_j^exp is the experimental value of coordinate j, and set their experimental values to the specified coordinates.

The comparison between the knee angle analyzed by Cortex data acquisition software and the knee angle analyzed by OpenSim musculoskeletal simulation software is shown in Fig. 7. The data in the figures represents the gait information of 1.5 cycles. The overall trend of knee joint angle obtained by two methods are similar and can be seen from the figures. However, there are still some differences. It can be seen that the variation range of knee joint angle obtained by using musculoskeletal simulation software OpenSim is larger and the variation trend is more stable. The reason of this phenomenon maybe is the Cortex software get the angles just by simple calculating with the collected position data. OpenSim combines the characteristics of musculoskeletal model in the calculation of joint angles, so the joint angles obtained by OpenSim are better than Cortex.

Inverse dynamics (ID)

The ID problem refers to: given the position q, velocity $\dot{q}$ and acceleration q of each joint of the robot at a certain moment, calculate the driving force (include: states or motion) imposed on each joint at this time. The ID can be solved by the Newton-Euler equation or the Lagrange equation.

Dynamics is the study of motion, the forces and torques that cause motion. The purpose of ID is to estimate the forces and torques required to produce a particular motion, and the results also used to predict how muscles contribute to motion. In order to calculate the forces and moments, the system’s equations of motion need to be solved iteratively. The motion equations are derived from the motion description and the mass property of the model. In the solution process, the IK is employed to calculate the joint angle and the ground reaction force during the experiment. And combined with the dynamic equilibrium conditions and boundary conditions, the net forces and moments at each joint are obtained.

When solving the ID problem, the data of motion and force measuring platform are needed to ensure that the number of equations of motion more than unknowns (degrees of freedom), which turns the problem into a statically indeterminate problem. The error of experimental motion data and the inaccuracy of the model will lead to Newton’s second law F=m*a invalid. In order to solve this dynamic discontinuity problem, residual forces and torques are introduced into a specific section of the model. The following equation is established, which relates the ground reaction moment to the residual moment. Where, ${\vec{F}}_{exp}$ is the ground reaction moment and ${\vec{F}}_{residual}$ is the residual moment. (2)${\displaystyle {\vec{F}}_{exp}+{\vec{F}}_{residual}=m\cdot \vec{a}}$

Inverse_Dynamics.sto: is generated by inverse kinetic operations, including time series, net joint moments of each bone joint, etc.

The net joints acting torque of the hip joints in three motion states are shown in Fig. 8. The net joint acting moment of hip joint is same to the angle of the hip joint, including adduction, flexion and rotation. Torque is drives of the body segments, therefore corresponds to the joint angles of the hip joint. The net joint torque of the knee joint is the joint torque connecting the thigh and the calf, which is significant in the research of exoskeleton rehabilitation robots of lower limbs.

Dynamics

Fourier is a rigid body robot, in the absence of friction and other disturbances, the dynamics of a serial n-link rigid robot can be written as: (3)${\displaystyle M\left(q\right)\ddot{q}+D\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=\tau +{\tau }_{OpenSim}.}$

We simplify Fourier to a 2-link rigid robot, where, q is the 2 × 1 vector of joint angle, τ is the 2 × 1 vector of joint torques, τ_OpenSim is the 2 × 1 vector of joint torques from the OpenSim software, M(q) is the 2 × 2 symmetric positive definite manipulator inertia matrix, $D\left(q,\dot{q}\right)$ is the 2 × 1 coriolis force and centrifugal force matrix, G(q) is the 2 × 1 gravity matrix. In this paper, τ_OpenSim as a reference torque input to the controller, to realize the anthropomorphic torque output of the exoskeleton robot.

The PD control law is given as flowing: (4)${\displaystyle \tau =K_d\dot{e}+K_{p}e-{\tau }_{OpenSim}}$

where, K_d and K_p are the 2 × 2 symmetric positive definite matrices, e = q_d − q, q_d is the desired joint angle. The τ_OpenSim is a bounded matrix.

${\displaystyle T_{\min }\le ∥{\tau }_{OpenSim}∥\le T_{\max }.}$ Eqs. (3) and (4) imply. (5)${\displaystyle M\left(q\right)\left({\ddot{q}}_d-\ddot{q}\right)+D\left(q,\dot{q}\right)\left({\dot{q}}_d-\dot{q}\right)+K_d\dot{e}+K_{p}e=0.}$

Then the $M\left(q\right)\ddot{e}+D\left(q,\dot{q}\right)\dot{e}+K_{p}e=-K_d\dot{e}$ is obtained. Considering the candidate Lyapunov function, (6)${\displaystyle V=\frac{1}{2}\dot{e^T}M\left(q\right)\dot{e}+\frac{1}{2}\dot{e^T}{K}_p\dot{e}}$

Eq. (6) is position definite, and the time derivative of the function, (7)${\displaystyle \dot{V}=\dot{e^T}M\ddot{e}+\frac{1}{2}\dot{e^T}\dot{M}\left(q\right)\dot{e}+\dot{e^T}{K}_{p}e.}$

There is an oblique symmetry feature: $\dot{M}\left(q\right)\text{-2}D\left(q,\dot{q}\right)=0$. Substituting the condition into Eq. (5), get (8)${\displaystyle \dot{V}=\dot{e^T}\left(M\ddot{e}+D\dot{e}+K_{p}e\right)=-\dot{e^T}{K}_d\dot{e}\le 0.}$

Then, the stability of design control system can be guaranteed.

Results

The simulation is carried out in MATLAB-Simulink, and the results are shown in Figs. 9 and 10. Figures 9A, 9C and 10A are the joints angle tracking and tracking error of OpenSim combined with PD controller, respectively. Figures 9B, 9D and 10B are the joints angle tracking and tracking error of the PD controller, respectively. Comparing Figs. 9A and 9B, it can be found that the PD controller combined with OpenSim has a better tracking effect, after the tracking, error-free tracking can be achieved, and the controller that only uses the PD control can clearly found that there is still an error in the peak position of the gait angle in the later stage of the tracking. In the method proposed in this paper, the joint torque from OpenSim plays a good role in compensating for the control of the controller. Comparing Figs. 10A and 10B, we can be found that the PD controller combined with OpenSim has a better tracking effect, the tracking error is smoother with only small fluctuations, and the stability is higher. The patient can enter the rehabilitation training comfortably, and achieve the safety and comfort of the rehabilitation training.

The parameter of PD controller, such as: Kp, Kd. Kp = diag(150, 150), Kd = diag(120, 120). The performance indexes of the two controllers designed are shown in Table 3.

The RMSE, ISE, and ITSE indicated that the two controllers are almost the same, suggesting that they have faster response speed and smaller oscillation. In terms of the properties of IAE and ITAE, the controller with OpenSim output torque is slightly better than the PD controller’s transient response and the transient response oscillation is smaller.