Dynamic Model of a Humanoid Exoskeleton of a Lower Limb with Hydraulic Actuators

Exoskeletons are the mechanical systems whose operation is carried out in close cooperation with the human body. In this paper, the authors describe a mathematical model of the hydraulic exoskeleton of a lower limb. The coordinates of characteristic points of the exoskeleton in the sagittal plane as a function of user height are presented. The mathematical models, kinematics, and kinetics equations were determined. The masses of the actuators and their dimensions were selected based on catalog data. The force distribution in the wearable system during the squat is shown. The proposed models allowed us to determine the trajectory of individual points of the exoskeleton and to determine the forces in hydraulic cylinders that are necessary to perform a specific displacement. The simulation results show that the joint moments depend linearly on actuator forces. The dynamics equations of the wearable system are non-linear. The inertia of the system depends on the junction variables and it proves that there are dynamic couplings between the individual axes of the exoskeleton.


Introduction
The development of robotics and computing power allows the design of devices for amplifying human power. This subject area is a present-day problem when we observe the processes of socio-cultural changes and the issues of society aging. In the 1950s, only 4.9% of the population lived more than 65 years and nowadays this percentage is about 20%, it is anticipated to be more than 35 % in the year 2050. During human aging, we become less able in terms of the body. Elderly people still want to be active, walk in the mountains, run, etc. The advantages of using modern devices can be seen in many areas of human performance. Exoskeletons are applied in medicine and rehabilitation by physiotherapists and patients. Some replicable movements with high precision can be made by machines which are applied in surgery or operating rooms. Mobile anthropomorphic robots are examples of the abovementioned machines which assist in the operation of human muscles and are called exoskeletons [1,2].
The first designs of devices for assisting human muscles were done in the 1950s at Argonne National Laboratory (USA) [3,4]. In 1965, the General Electric company started investigations on the exoskeleton Hardiman I (Human Augmentation Research and Development Investigation Man) [5,6]. The mass of this device was equal to 680 kg and one upper limb could raise a load of 340 kg. At the end of the 1980s, Jeffrey Moore from Los Alamos Laboratory described a vision of a combination of an exoskeleton with a spacesuit [7].
During the last 20 years, exoskeletons have been developed considerably [8]. This could be done due to the grant aid from military sources like the Defense Advanced

Exoskeleton Scheme
The model of a hydraulic exoskeleton with three actuators is shown in Figure 1a. A kinematic description and mathematical model with the movement simulation are presented in [37].

Coordinates
The origin of the coordinate system was placed in the ankle joint to determine the forces for rising of the human body from the squat position to the standing position. A kinematic model is shown in Figure 1b. Transformation matrixes can be written in the following form: Locations of characteristic points can be expressed by the following dependence: where SC l -location of center of gravity for shank, SC t -location of center of gravity for thigh, J † -location of center of gravity for third actuator, P † -location of center of gravity for first actuator. If the following parameters are known: angle of initial position (squat), angle of final position (standing position), velocity and acceleration values for start and end of the movement, and time of movement. Then, we can determine the locations of individual points with the application of the forward kinematics [37]. The sum of individual angles in the hip joint θ h , knee joint θ k , and ankle joint θ a and the sum of individual velocitiesθ h ,θ k ,θ a and angular accelerationsθ h ,θ k ,θ a can be defined by If c and s denote cosine and sine functions then the coordinates of mass centers, velocities, and accelerations of the single elements of an exoskeleton with a lower limb can be described by the equations The movement trajectories of mass centers of the trunk, thigh, shank, and actuators 1, 2, and 3 are shown in Figure 2. Equation (4) was applied in order to determine the trajectory, where t 0 -initial position (squat), t k -standing position. The arrows indicate the movement direction of mass centers in relation to the coordinate system origin. The highest value of displacement on the x-axis can be observed for the trunk and its value is −0.28 m. The highest value of displacement on the y-axis occurs for actuator 3 which drives the ankle joint and its value is 0.2 m. The lowest value of displacement on the x-axis can be observed for the mass center of the shank, its value is 0.072 m-this fact is logical because the distance from the origin of the coordinate system is small.  The velocity values on x-and y-axis for individual point masses are presented in Figure 3. The highest values of velocity (about 1 m/s) were obtained by the actuator for the knee joint on the y-axis. The lowest values of velocity on the x-and y-axis were observed for the mass center of the shank.

Time [s]
Time Differentiation of Equation (5) gives us the following acceleration components: The values of acceleration components of mass centers for individual elements of the model are presented in Figure 4.

Time [s]
Time [s] The component values of angles, velocities, and angular accelerations for single joints during a change in position from squatting to standing upright in 1s are shown in Figure 5a. The displacement, velocity, and acceleration of actuators that drive the joints are illustrated in Figure 5b.  Force distribution in hydraulic exoskeleton during knee bending is shown in Figure 1c. The single parts of limbs and the exoskeleton (these are rigid bodies with homogenous mass distribution) are replaced by the finite set of point masses which are aggregated in the pre-selected points [38]. To determine the characteristic angles, we should use the geometrical dependencies (7) and pay attention to the selection of proper constant distances between the characteristic points of the limb and exoskeleton. According to Jezierski, the angle for this type of drive is 180 • [39]. For angle values of 0 • and 180 • , a problem exists because, for these values, we cannot drive the joint. To prevent this situation, we should use a structure called the four-bar linkage mechanism and we should take into account the predicted load of mechanisms and the selection of proper parameters [40].
where the value of angles γ 1 , γ 2 is constant, whereas α 1 , α 2 and the distance between points BP, EI are functions of angles in hip and knee joints θ h , θ k .

Results
The diagram of forces, moments of inertia, and reactive forces for free limbs is presented in Figure 6.  Inclination angles of straight lines of actuator operation ξ F1 , ξ F2 , ξ F3 są F 1 , F 2 , F 3 can be calculated by the following equation (we should apply here the coordinates from Equation (2)): To determine the distances between the individual masses with characteristic points (hip, knee, and ankle joint) in the function of angles θ h , θ k , θ a , we should use the dependencies (4) and (2), and then we get the vectors: where the absolute value of a vector with characteristic points G, K can be calculated as a root of the sum of the component squares x and y: The other values like r (bC) , r (tG) , and r (lK) are constant and depend on the user's height. The components of force F 1 , F 2 , and F 3 are projected on the individual axes of a pre-selected coordinate system and have the following form: Driving moments M F1 (F 1 ,α 1 ), M F2 (F 2 ,α 2 ), and M F3 (F 3 ,α 3 ) can be calculated as where the forces F 1 , F 2 , and F 3 from actuators 1, 2, and 3 can be calculated by The moments M F1 (θ h ), M F2 (θ k ), and M F3 (θ a ) depend in a linear way on forces F 1 , F 2 , and F 3 which are generated by the actuators, and depend in nonlinear way on angles α 1 , α 2 , and α 3 . The equilibrium equation of the system presented in Figure 6 in a given time can be written in the following way (after some transformations): where we can use the equation of state to determine the forces F 1 , F 2 , and F 3 , reactions C, G, and K, and moments in characteristic points of the system: Mass, dimensions, and location of the center of gravity for actuators have been defined with information from catalogs [41][42][43]. Every actuator was modeled in SolidWorks and this allowed us to calculate the moments of inertia [44]. The dependencies (15) are valid for any mass and length of a limb part with the individual elements of the exoskeleton. The orientation of the coordinate system is given in Figure 6c. The movement trajectory of the exoskeleton was divided into one hundred elements where the location of characteristic points, components of forces, and moments were defined. We define the equilibrium conditions by using Equation (15) and next we calculate the forces in every actuator. From the first equation, we can calculate moment M Gb with point C-it results from the displacement of the center of gravity of the upper part of the body with an exoskeleton in relation to the hip joint. Next, we can determine the force F 1 and the reaction of the hip joint C x and C y which are necessary to balance the moment. The next step is the determination of the moment in the knee joint G (we need to take into consideration the moment from the hip joint), force F 2 for the equilibrium of the system, and reactions G x and G y . The last step is the calculation of the moment in the ankle joint K, force F 3 in the actuator, and reactions K x , K y and M x , M y . The abovementioned procedure is repeated for every new position. Force values F 1 , F 2 , and F 3 for actuators with components x and y and reactions in individual joints are presented in Figure 7. The force value F 1y for the final position is 9 N. It results from the final position of the inclination angle of the trunk which is 16 • , the center of gravity of the trunk is displaced forwards about 2.8 cm with the hip joint. Component reaction G x in the hip joint for angle range θ h = 85 • − 45 • has a negative value because of the action line of force F 1 . Longitudinal axes of actuators are not parallel to individual limb elements and the rotation of the limb causes changes in the location of mass centers in relation to individual points. Human muscles do not give any forces in the analyzed model. In reality, the forces in flexor muscles and extensors can considerably reduce the participation of actuators in keeping the system equilibrium. We should note that the force in actuator 2 is acting both on point E which is connected with the thigh part of the exoskeleton and point I which is connected with the shank part.

Time [s] Time [s] Time [s]
Moments of inertia of human body parts in relation to the center of gravity were calculated based on the anthropometric data-the mass of the exoskeleton was divided and assigned to every limb (i.e., trunk with the exoskeleton, head, upper limbs). For a man with a mass of 80 kg and a height of 175 cm, the values of inertia moments are the following: I b = 5.6257 kgm 2 , I t = 0.166 kgm 2 , and I l = 0.0787 kgm 2 . Moments of inertia for actuators were calculated for 3D models in SolidWorks [44]. Based on the movement trajectory of individual elements, we calculated the values of velocity and acceleration for every point in the system.
The equations of the translational motion of a hydraulic exoskeleton with a limb in the sagittal plane and the equations of rotational motion in relation to hip, knee, and ankle joints can be described by using dependencies (6) in the following way: where moments of inertia depend on the distance between point masses (Equation (4)) and the joint. These can be calculated as: Dynamic equations of motion (16) can be described in a general matrix form as where A N -matrix of drive mass, A B -matrix of mass of human body parts, I-matrix of moments of inertia. The individual equations were transformed by using dependencies (12)- (15) due to the size of each abovementioned matrix. On the left side of nine equations, we marked forces and moments of forces as unknown values: where r BCx,y , r PGx,y , r CGx,y , r JKx,y are components x and y for the distances between points BC, PG, CG, JK. The forces in point M (fixing of piston rod of cylinder 3) are M x = F 3x and M y = F 3y . The procedure is the same as in the case of the analysis of static equilibrium of the system. When writing the equations in a computer program, we should pay particular attention to sign correctness where the positive moment is assumed for a right-handed screw rotation. In the dynamic model, we have to take into account moments of inertia that are connected with mass distribution in the body. The simulation results are illustrated in Figure 8. We can observe that the values of forces and reactions are the same for dynamic and static models-during the start and stop of movement for characteristic points.
The system inertia and time of displacement from the initial position to the final one have a deciding impact on necessary forces for acceleration and deceleration of the system. In the dynamic model, the force from human muscles is zero. As was mentioned earlier, the dynamic equations of the hydraulic exoskeleton of a lower limb are nonlinear. The system inertia depends on junction variables and it proves the occurrence of dynamic coupling between individual axes of the exoskeleton. The essential role here is played by the velocity coupling, i.e., centrifugal force and Coriolis force. The current configuration of the exoskeleton influences the gravitational force which affects the device dynamics.

Discussion
In this work, the authors developed a hydraulic model of an exoskeleton lower limb. We built a geometrical model of a hydraulic exoskeleton based on anthropometrical parameters. The paper presents the mathematical models related to the dynamics of the exoskeleton based on the matrix notation. The proposed mathematical models are built as a function of human height. This method allows for determining the location of characteristic points of the wearable system. This is especially important when designing this type of device. The starting and ending position of the wearable device was calculated with the determination of the movement trajectory based on the method presented in [37]. The proposed mathematical model of each part of the exoskeleton takes into account the weight of an element of the wearable device and each part of the user's body. The moments of inertia of the human body were based on the literature data. Static analysis of exoskeletons was done in every position, from the initial position (squat) to the final position (standing position). Next, the time of movement was determined. The necessary forces and moments which are needed to make the displacement were determined on the basis of the presented dynamics equations and including the moments of inertia for every element. The correctness of the models was verified by comparing the consistency of the force values at the beginning and the end of the motion with the static and dynamic analysis.
Of course, the proposed model has some limitations. One is the fact that it is developed in the sagittal plane. The real model has not been developed yet, therefore, it will be necessary to make a prototype and verify the proposed models. Much work needs to be done to accurately design the detailed exoskeleton model, analyze its performance, and make a real model. In future work, the authors will collect data from the real wearable system which will be built. To optimize the operation of the exoskeleton, it is also necessary to conduct parametric studies to select the best possible values of the design variables. The current results did not consider the change in actuator mass when filling the actuator chamber. θ t thigh angle to coronal plane, raḋ θ h hip joint angular velocity, rad/ṡ θ k knee joint angular velocity, rad/ṡ θ a ankle joint angular velocity, rad/s θ h hip joint angular acceleration, rad/s 2 θ k knee joint angular acceleration, rad/s 2 θ a ankle joint angular acceleration, rad/s 2 M h , M k , M a cosine matrices, Nm A N matrix of drive mass, kg A B matrix of mass of human body parts, kg I matrix of moments of inertia, kgm 2