Optimization Based Dynamic Human Motion Prediction with Modular Exoskeleton Robots as Interactive Forces: The Case of Weight Lifting Motion

Abstract

The optimization-based dynamics model is formulated for the weight lifting motion with human and exoskeleton model as interactive force term in this chapter. In the optimization algorithm, the human motion is defined as variables so that the motion which we want to generate (box lifting motion in this case) can be predicted. The objective function or cost function is defined as performance measure which can be switched by developer. In this paper we use the summation of each joint torque square which is considered as the dynamic effort for the motion. Constraints are defined as joint limits, torque limits, hand position, dynamic balance, exoskeleton assistive points, etc. Interaction force form exoskeleton robot can be derived as generalized coordinates and generalized force which are related to inertial reference frame and human body frame. The results can show how effective the exoskeleton robots are according to their assistive force.

1. Introduction

To design or to assess the exoskeleton robot, it is necessary to simulate human motion with interactive force from exoskeleton robot. Thus, we need the human modeling and simulation which method the interactive force from exoskeleton can be applied at. This is key motivation in this study and we formulate human motion simulation with modular exoskeleton model and predict human motion with interactive from modular waist and knee exoskeleton robot. To do that we formulate human motion with assistive force from exoskeleton robot. There are many good researches based on optimization techniques for simulation in different areas [1, 2, 3, 4]. There are couple of human modeling and simulation software such as OpenSim [5]. By given motion OpenSim can generate muscle forces for that motion. Other is Santos which is optimization-based motion simulation so called predictive dynamics [6]. Some of human motion simulation is developed under optimization-based motion simulation [7, 8]. And, lifting motions are studied with optimization based technique as well [9, 10]. The modeling and simulation for human-exoskeleton is developed in [11]. In this study, we also use optimization-based motion prediction and simulation method with recursive Lagrange’s equations of motion which is highly nonlinear. It provides predicted human motion as well as joint torques and ground reaction forces for weight lifting. This method also gives us pretty fast calculation time and accurate results.

2. Mechanical modeling

2.1 Modular exoskeleton robot

There are many researches about exoskeleton with human-robot cooperation [12]. Due to the wearability, convenience, comfort and easy portability, modular exoskeleton robot is becoming a trend in the industrial work environment nowadays such as construction site, heavy industry, medical care, logistics, maintenance, manufacturing process, etc. In this study, the modular means the modular type according to the body parts of human being. For example, shoulder modular exoskeleton, knee modular exoskeleton. Biggest merit of modular exoskeleton may be the bringing more comfortability rather than full-body exoskeleton robots. Of course, it is possible only in the industrial area. If we look for the purpose of rehabilitation in hospital, it may be different story. Also, once we narrow down the application area, modular exoskeleton robot can be lighter, have more simple structure and can be more compact. Some area, you do not need active exoskeleton robot, and just passive exoskeleton is fine. Also, modular exoskeleton robot is applicable either together or separate case by case. Furthermore, modular exoskeleton is more economical compare to the full-body exoskeleton robot. The following Figure 1 shows the concept design of modular exoskeleton robot which we are currently developing. In this study, exoskeleton assistive force can be applied human body as an external force through the optimization process.

Figure 1.

The concept design of modular exoskeleton robot.

2.2 Human model

Human model is constructed using mechanical structure with revolute joints and rigid body. Thus, each body segment is connected with revolute joint. There are 49 degrees of freedoms (DOF) in each body joint. Such as knee is 1 DOF, shoulder is 3DOFs etc. There are 6 DOFs for global translation and rotation for human system to the inertial reference frame. We put inertial reference frame at the point between the foot on the ground as human model is standing. We assume that body segment is rigid body and there is no muscle and tissue in this study. Therefore, it is assumed that all muscle force is converted to the joint torque. Also, we used GEBOD software to generate dynamic properties of body segment for example thigh, pelvis, and torso [13]. In this study we used 50 percentile male data which is representing average male size. Figure 2 describes the mechanical structure of current human model to generate weight lifting motion with waist and knee modular exoskeleton robots. Each z-axis has DOF and transformation matrix is combined sequentially from inertial reference frame to the head, hands, and toes as a branch. The virtual branch depicts global DOFs which is mentioned in previous – 3 global translations and 3 global rotations. The torso part of the human model has 4 spine joints and there are total 12 DOFs. Then, it leads to the right arm branch, left arm branch, and head branch. Right arm and left arm branch has 4 joints and 9 DOFs respectively including clavicle joint. Head branch has 2 joints and 5 DOFs. Right leg and left leg has 4 joints and 7 DOFs respectively.

Figure 2.

Mechanical structure of human model.

3. Kinematics and dynamics analysis

3.1 Kinematics

Denavit-Hartenberg (DH) method is used to analyze the kinematics of human motion [14]. In the Denavit-Hartenberg method any point ⁱr can be transferred to the global reference frame as ⁰r in Figure 3 and it can be presented as Eq. (1).

Figure 3.

Articulated chain.

$${}^0\mathbf{r}={{}^0\mathbf{T}}_i{}^i\mathbf{r}$$E1

Transformation matrix ⁰T_i in Eq. (1) can be obtained as follows.

$${{}^0\mathbf{T}}_i={{}^0\mathbf{T}}_1{{}^1\mathbf{T}}_2\cdots {{}^{i-1}\mathbf{T}}_i=\prod _{n=1}^i{{}^{n-1}\mathbf{T}}_n$$E2

$${{}^{i-1}\mathbf{T}}_i=\left[\begin{array}{cccc}cos{q}_i& -cos{\alpha }_{i}sin{q}_i& sin{\alpha }_{i}sin{q}_i& a_{i}cos{q}_i\\ sin{q}_i& cos{\alpha }_{i}cos{q}_i& -sin{\alpha }_{i}cos{q}_i& a_{i}sin{q}_i\\ 0& sin{\alpha }_i& cos{\alpha }_i& d_i\\ 0& 0& 0& 1\end{array}\right]$$E3

where the parameters q, α, a, d, are DH parameters. The q_i is the joint angle between the x_i-1 axis and the xi axis about the z_i-1 axis according to the right hand rule. The d_i is the distance between the origin of the i-1th coordinate frame and the intersection of the z_i-1 axis with the x_i axis along the z_i-1 axis. a_i is the distance between the intersection of the z_i-1 axis with the x_i axis and the origin of the i^th frame along the x_i axis. α_i is the angle between the z_i-1 axis and the z_i axis about the x_i axis according to the right hand rule. In here, we use q_i as our generalized coordinates. As mentioned before, the inertial reference frame is located at the point between foot. The origin of inertial reference frame is O in the above Figure 2. Thus, all kinematic chain starts from origin of inertial reference frame. For the efficiency of calculation time, we use recursive way for kinematic information as follows:

$${\mathrm{A}}_i={\mathrm{A}}_{i-1}{\mathrm{T}}_i$$E4

$${\mathrm{B}}_i={\dot{\mathrm{A}}}_i={\mathrm{B}}_{i-1}{\mathrm{T}}_i+{\mathrm{A}}_{i-1}\frac{\partial {\mathrm{T}}_i}{\partial q_i}{\dot{q}}_i$$E5

$${\mathrm{C}}_i={\dot{\mathrm{B}}}_i={\mathrm{C}}_{i-1}{\mathrm{T}}_i+2{\mathrm{B}}_{i-1}\frac{\partial {\mathrm{T}}_i}{\partial q_i}{\dot{q}}_i+{\mathrm{A}}_{i-1}\frac{{\partial }^2{\mathrm{T}}_i}{\partial q_i^2}{\dot{q}}_i^2+{\mathrm{A}}_{i-1}\frac{\partial {\mathrm{T}}_i}{\partial q_i}{\ddot{q}}_i$$E6

where A_i is position matrix, B_i is velocity matrix, and C_i is acceleration matrix.

3.2 Dynamics

Once we obtain the kinematic information, we can use them to calculate dynamics of motion simulation. The equations of motion are derived from Lagrange’s equation.

$$Q_i=\frac{d}{\mathit{dt}}\left(\frac{\partial L}{\partial {\dot{q}}_i}\right)-\frac{\partial L}{\partial q_i}$$E7

where L is Lagrangian which is L = K – V (kinetic energy – potential energy), and q_i is the generalized coordinate. Total kinetic energy is derived as

$$K=\sum _{j=1}^n{K}_j=\frac{1}{2}\sum _{i=1}^n\mathit{tr}\left({\dot{\mathbf{A}}}_i{\mathbf{J}}_i{\dot{\mathbf{A}}}_i^T\right).$$E8

where A_i is joint angle matrix and J_i is inertia matrix at i^th reference frame. Also, the potential energy for system can be given as

$$V=-\sum _{j=1}^n{m}_j{\mathbf{g}}^T{\mathbf{A}}_j{{}^j\overline{\mathbf{r}}}_j-\sum _{j=1}^n{\mathbf{f}}_k^T{\mathbf{A}}_j{{}^k\mathbf{r}}_f{\delta }_{\mathit{jk}}$$E9

where m_j is the mass of the body segment represented as j^th reference frame, g is gravity force vector, f_k is an external force which is defined in global reference frame and acting on the body segment expressed in k^th reference frame, ^kr_f is the location of external force acting on the link expressed in the k^th reference frame, δ_jk is Kronecker delta. Then, the equations of motion can be derived from above equations (Eqs. (7)–(9))

$${\tau }_i=\mathit{tr}\left(\frac{\partial {\mathbf{A}}_i}{\partial q_i}\sum _{j=i}^n{{}^i\mathbf{T}}_j{\mathbf{J}}_j{\ddot{\mathbf{A}}}_j^T\right)-{\mathbf{g}}^T\frac{\partial {\mathbf{A}}_i}{\partial q_i}\sum _{j=i}^n{m}_j{{}^i\mathbf{T}}_j{{}^j\overline{\mathbf{r}}}_j-{\mathbf{f}}_k^T\frac{\partial {\mathbf{A}}_i}{\partial q_i}\sum _{j=i}^n{{}^i\mathbf{T}}_j{{}^{k}r}_f{\delta }_{\mathit{jk}}-{\mathbf{G}}_i{\mathbf{A}}_{i-1}{\mathbf{z}}_0$$E10

where τ_i is the joint torque acting on the joint represented with generalized coordinate q_i, G_i is external moment which is defined in the global reference frame and z₀ is [0 0 1 0]^T. Ground reaction force can be calculated using global force transformation the mechanical system of human body. Obtained joint torques are used to evaluate joint torque constraints and objective functions in later section.

4. Optimization formulation

Optimization based motion simulation in this chapter is performed according to the following process:

Step1: Prepare input data.

Step2: Function approximations for joint variables.

Step3: Kinematics analysis.

Step4: Dynamics analysis.

Step5: Objective function evaluation.

Step6: Constraints evaluation.

Step7: Print out if converge. Otherwise go back to step 2.

Kinematics and dynamics are covered in previous section, so we will discuss the joint variables, objective function, and constraints evaluation in this section.

4.1 Variables

To generate weight lifting motion, our optimization variable is joint angle profiles. These joint angle profiles are approximated using B-spline function approximation [15]. The control points in B-spline function approximation are updated in each iteration through the optimization process. Also, we used clamped B-spline which is that the starting point and end point of spline curve are matching to the control points. Following equation describes the B-spline approximation for the joint angle profiles q_i.

$$q_i\left(\mathbf{t},\mathbf{p}\right)=\sum _{j=0}^m{N}_j\left(\mathbf{t}\right)p_j$$E11

where t = {t₁,…,t_s} is knot vector, p = {p₁,…,p_m} is the control points, N_j is the basis function of B-spline function approximation. Here, control points P becomes optimization variables.

4.2 Performance measure and objective function

We define the energy consuming (dynamic effort) for the given motion as the performance measure for weight lifting motion prediction and simulation. In this study, we use the joint torque square which is proportional to the mechanical energy. This mechanical energy is a reasonable criterion to minimize [4]. Then, it is formulated as objective function in optimization formulation as follows:

$$f={\int }_0^T\left[\sum w_1{{\tau }_i}^2\right]\mathit{dt}$$E12

where w_i is the weighting parameters of each joint, τ_i is the joint torque of each joint. The joint torque is calculated from above dynamics equation (Eq. (10)).

4.3 Constraints

Constraints are the one of the motion control way in this optimization formulation and we used minimal set of constraints to generate motion in this formulation. The list of constraints are as follows:

1. Joint angle limits

2. Joint torque limits

3. Foot position and hand position

4. Stability condition

The joint angle limits and torque limits for the human motion are determined based on the literature [16, 17, 18, 19]. Zero Moment Point(ZMP) method is used for the stability condition constraint [20]. Each constraint is formulated as follows accordingly:

$${\mathbf{q}}^L\le \mathbf{q}\left(t\right)\le {\mathbf{q}}^U$$E13

$${\mathbf{\tau }}^L\le \mathbf{\tau }\left(t\right)\le {\mathbf{\tau }}^U$$E14

$${\mathbf{r}}_i\left(t\right)={\tilde{\mathbf{r}}}_i$$E15

$$z_{\mathit{zmp}}\in \mathrm{FSR},x_{\mathit{zmp}}\in \mathrm{FSR}$$E16

4.4 Analytical gradients

The analytical gradients of each constraints and objective function are provided to the optimization solver. These analytical gradients improve the accuracy as well as calculation time for convergence in optimization process. Analytical gradients for constraints and objective function are as follows:

$$\begin{array}{l}\frac{\partial C}{\partial x_1}=\frac{\partial C}{\partial q_1}\frac{\partial q_1}{\partial x_1}+\frac{\partial C}{\partial q_2}\frac{\partial q_2}{\partial x_1}+\cdots +\frac{\partial C}{\partial q_n}\frac{\partial q_n}{\partial x_1}+\frac{\partial C}{\partial {\dot{q}}_1}\frac{\partial {\dot{q}}_1}{\partial x_1}+\cdots +\frac{\partial C}{\partial {\dot{q}}_n}\frac{\partial {\dot{q}}_n}{\partial x_1}+\cdots +\frac{\partial C}{\partial {\ddot{q}}_n}\frac{\partial {\ddot{q}}_n}{\partial x_1}\\ \frac{\partial C}{\partial x_2}=\frac{\partial C}{\partial q_1}\frac{\partial q_1}{\partial x_2}+\frac{\partial C}{\partial q_2}\frac{\partial q_2}{\partial x_2}+\cdots +\frac{\partial C}{\partial q_n}\frac{\partial q_n}{\partial x_2}+\frac{\partial C}{\partial {\dot{q}}_1}\frac{\partial {\dot{q}}_1}{\partial x_2}+\cdots +\frac{\partial C}{\partial {\dot{q}}_n}\frac{\partial {\dot{q}}_n}{\partial x_2}+\cdots +\frac{\partial C}{\partial {\ddot{q}}_n}\frac{\partial {\ddot{q}}_n}{\partial x_2}\\ ⋮\\ \frac{\partial C}{\partial x_m}=\frac{\partial C}{\partial q_1}\frac{\partial q_1}{\partial x_m}+\frac{\partial C}{\partial q_2}\frac{\partial q_2}{\partial x_m}+\cdots +\frac{\partial C}{\partial q_n}\frac{\partial q_n}{\partial x_m}+\frac{\partial C}{\partial {\dot{q}}_1}\frac{\partial {\dot{q}}_1}{\partial x_m}+\cdots +\frac{\partial C}{\partial {\dot{q}}_n}\frac{\partial {\dot{q}}_n}{\partial x_m}+\cdots +\frac{\partial C}{\partial {\ddot{q}}_n}\frac{\partial {\ddot{q}}_n}{\partial x_m}\end{array}$$E17

$$\frac{\partial \tau }{\partial x_k}=\frac{\partial \tau }{\partial q_j}\frac{\partial q_j}{\partial x_k}+\frac{\partial \tau }{\partial {\dot{q}}_j}\frac{\partial {\dot{q}}_j}{\partial x_k}+\frac{\partial \tau }{\partial {\ddot{q}}_j}\frac{\partial {\ddot{q}}_j}{\partial x_k}$$E18

where C represents constraints, q_i is joint angle profile, x_i is control point of B-spline function approximation, τ is joint torque profile. The analytical gradients of each constraint are calculated in the form of Eq. (17). For the analytical gradient of objective function is calculated by using Eqs. (12) and (18).

5. Numerical experiments

5.1 Weight lifting motion

Figure 4 depicts weight lifting motion. The weight is defined as W and it is applied as and external load in dynamics equilibrium equations of motion. Weight object is virtual and only hand location is guided by position constraints. Foot is located on the ground and weight is moving up vertically in 0.15 m from 0.5 m above ground and away from heal position by 0.25 m as shown in Figure 4. Dynamic stability constraint which is ZMP are imposed on the foot support region area. Joint angle limits and joint torque limits are imposed based on the literature review. It is assumed that the assistive moment is applied hip joint and knee joint which axes are parallel to the horizontal axis of inertial reference frame.

Figure 4.

Weight lifting motion with waist and knee exoskeleton robots.

5.2 Simulation setting

For the optimization, sequential quadratic programming (SQP) is adopted. The sequential quadratic programming is very effective for the large scale nonlinear constrained optimization problem. The commercial software SNOPT is used which is well known as the effectiveness for large scale nonlinear problem SQP solver [21]. Total number of variables for optimization is 330 and total number of constrains is 1,766. The weighting parameters in objective function are set to equally in current study. Numerical experiment is performed while the weights are set to 15 kg and 30 kg. Also, different assistive moments from exoskeleton robot are tested in the experiment. The postprocessing for animation and snapshot was performed using Commercial software MATLAB.

6. Results

The simulation results are shown in following Figure 5. The assistive force from exoskeleton robot is tested from 10 N/m to 50 N/m for waist and knee exoskeleton. Then, we checked energy consumption for the lifting motion from the simulation. The results are Table 1 and Figure 6. As shown in the Table 1, it is obvious that the exoskeleton robot reduces the total energy consumption of weight lifting motion. Most energy minimized case for each 15 kg and 30 kg lifting case is written in italic in the table. In some case, the total energy is more than no assistive force case. It might be the assistive force is bothering the balance mechanism of human body while weight lifting motion so human uses more energy to recover the balance back.

Figure 5.

Simulation snapshot for weight lifting motion.

Lifting weight	Exo assistive force (N/m)		Objective function value (energy consumption)	Lifting weight	Exo assistive force (N/m)		Objective function value (energy consumption)
	waist	knee			waist	knee
15 kg	0	0	3.3260043054–002	30 kg	0	0	6.6957529363–002
	10	10	3.0190845608–002		10	10	4.6085529919–002
	20	10	2.7882680491–002		20	10	5.8372796346–002
	30	10	2.6598938935–002		30	10	5.3928838227–002
	40	10	2.6341574623–002		40	10	5.0509039531–002
	50	10	2.7108828860–002		50	10	4.8111291575–002
	10	20	2.4055151982–002		10	20	4.5952969619–002
	20	20	2.7788637489–002		20	20	4.5957540847–002
	30	20	2.6149510912–002		30	20	5.4008290188–002
	40	20	2.5534968466–002		40	20	5.0232961805–002
	50	20	2.5946755651–002		50	20	4.7482075385–002
	10	30	2.3808412818–002		10	30	4.6126512561–002
	20	30	2.3808412818–002		20	30	4.3911156137–002
	30	30	2.6053418696–002		30	30	4.4983526018–002
	40	30	2.5083330659–002		40	30	5.0310773358–002
	50	30	2.5137987323–002		50	30	4.7204074709–002
	10	40	2.4208434896–002		10	40	12.775145795–002
	20	40	2.2956322073–002		20	40	4.6699616710–002
	30	40	2.2447356818–002		30	40	4.6284398370–002
	40	40	2.4985189230–002		40	40	4.2840986427–002
	50	40	2.4684139732–002		50	40	4.7280245854–002
	10	50	2.4132838565–002		10	50	16.305738815–002
	20	50	2.3243927837–002		20	50	10.469578598–002
	30	50	2.2435740271–002		30	50	4.7176051325–002
	40	50	2.1731970674–002		40	50	4.3894974761–002
	50	50	2.4583949090–002		50	50	4.2854122138–002

Table 1.

Numerical experiment results for weight lifting motion.

Figure 6.

The contour plot of numerical experiment results.

7. Conclusions

We have studied optimization-based motion simulation with modular waist and knee exoskeleton robot as an assistive force. We used Denavit-Hartenberg method for kinematics, Lagrange’s equations of motion with external force and moment term, B-spline function approximation. In motion simulation, the performance measure is mechanical energy which is presented as the summation of joint torque squares. Minimal constraints are applied such as joint angle limits, torque limits, dynamic balance, and hand/foot positions., the optimization process find out the minimized energy consumed motion under the assistive forces from the modular waist and knee exoskeleton robots which are applied during the weight lifting motion.

This method provides unique feature with human-exoskeleton modeling and simulation area. It can give us predictive motion of human so that the exoskeleton parameters are adjusted based on the predicted motion simulation. Also, human motion can be generated automatically for the control algorithm of robot to collaborate with human. It can be used as evaluation and assessment tool for the design parameters of exoskeleton robot development in any given tasks according to human factors. Furthermore, this can reduce the development cost of exoskeleton because not many prototypes are necessary and provides safe design and test process during the exoskeleton development procedure. Of course, the musculoskeletal model should be developed for more accurate calculation of human factors and it will be remained as future works.

Acknowledgments

This research was financially supported by the Institute of Civil Military Technology Cooperation funded by the Defense Acquisition Program Administration and Ministry of Trade, Industry and Energy of Korean government under grant No. 19-CM-GU-01.