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

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.


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.

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.

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.

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

Figure 2.
Mechanical structure of human model.
Transformation matrix 0 T i in Eq. (1) can be obtained as follows.
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: where A i is position matrix, B i is velocity matrix, and C i is acceleration matrix.

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.
where L is Lagrangian which is L=K-V (kinetic energypotential energy), and q i is the generalized coordinate. Total kinetic energy is derived as 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 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, k r 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)) 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 0 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.

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.
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.

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 .
where t ={t 1 , … ,t s } is knot vector, p ={p 1 , … ,p m } is the control points, N j is the basis function of B-spline function approximation. Here, control points P becomes optimization variables.

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: 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)).

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: 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:

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: where C represents constraints, q i is joint angle profile, x i is control point of Bspline 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). 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.

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.

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.

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.