Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons

i Equation Chapter 1 Section 1 End of Master Project Master of Industrial Engineering Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons Author: Guillermo Muñoz Pérez Tutors: Juana Mª Mayo Núñez Joaquín Ojeda Granja Mechanical Engineering and Manufacturing Department Higher Technical School of Engineering University of Seville Seville, 2021 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons iii End of Master Project Master of Industrial Engineering Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons Author: Guillermo Muñoz Pérez Tutores: Juana Mayo Núñez Joaquín Ojeda Granja Dpto. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons v End of Master Project: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons Autor: Guillermo Muñoz Pérez Tutores: Juana Mª Mayo Núñez Joaquin Ojeda Granja El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons vii Agradecimientos Tras tantos meses de duro trabajo, videoconferencias, horas frente al ordenador, lágrimas de agobio y desesperación, y momentos de alegría tras cada pequeño logro conseguido, este Trabajo llega a su fin, y con él mi etapa académica, que no formativa, pues uno siempre se está formando y aprendiendo en este camino que es la vida. Y mucho he aprendido en estos más de 20 años. Empezando por la señorita Antoñita (descanse en paz), de Educación Infantil, que me enseñó a contar y a escribir las letras y los números en un papel. Siguiendo por las valiosas lecciones de esos Maestros y Profesores que dignifican el nombre de su profesión, y no puedo evitar homenajear a mi recordado Maestro Don Fausto (descanse en paz). Allá donde esté, este Trabajo va dedicado a usted. Y concluyendo con los valiosos consejos que me han aportado mis Tutores en este Trabajo, Juana y Joaquín, de los que tanto he aprendido y que trataré siempre de poner en práctica en mi presente y futuro personal y profesional. Pero por supuesto, debo citar a mi familia. Esa familia que me ha criado, cuidado, educado y formado como persona. Esos padres, Frank y Maricarmen, que se han desvivido por mí y por mis hermanos, y esos hermanos, Javi y Carmen, que siempre están ahí para lo bueno y para lo malo. Y por último, pero no menos importante, sino todo lo contrario. La familia que he elegido, la persona que me ha acompañado desde el principio hasta el fin de este Trabajo, tragándose todos mis lamentos y agobios, apoyándome en todo momento, aportándome los consejos más valiosos del mundo y alegrándose de mis éxitos y avances tanto como de los suyos. Azahara, amor mío, este Trabajo, y la vida en general, serían mucho más difíciles sin ti. Te quiero. Guillermo Muñoz Pérez Sevilla, 2021 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons ix Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons xvii LIST OF TABLES Table 3.1. Experimental parameters in passive force-length equation by De Groote-Fregly [19] 30 Table 3.2. Parameters in force-velocity equation by De Groote-Fregly [19] 31 Table 5.1. Moment of inertia of each bar 50 Table 5.2. Main muscle properties. Data from OpenSim muscles [6] 54 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons xix LIST OF FIGURES Figure 1.1. Industrial worker wearing an Exhauss passive exoskeleton (Theurel et al., [3]) 1 Figure 1.2. active exoskeleton worn in lower limbs (Yale University, [4]) 1 Figure 1.3. Scheme of single shooting problem. [12] 2 Figure 2.1. Operating scheme of OpenSim Moco. [13] 4 Figure 2.2. Lineal interpolation through nodes in optimal control problem. Kelly, M. [15] 8 Figure 2.3. Operating scheme of MocoInverse problems 14 Figure 2.4. Operating scheme in MocoTrack problems 15 Figure 2.5. Operating scheme of MocoStudy problems 16 Figure 3.1. 3-element Hill’s model. Martínez Reina, J. [21] 19 Figure 3.2. Parameters in Hill’s muscle model. García Vallejo [25] 20 Figure 3.3. Activation and deactivation time constants. [25] 21 Figure 3.4. Tendon force length curve. OpenSim Documentation [28] 23 Figure 3.5. Passive force length curve. OpenSim Documentation [28] 24 Figure 3.6. Active force length curve. OpenSim Documentation [28] 24 Figure 3.7. Tendon force length curve. OpenSim Documentation [28] 25 Figure 3.8. Tendon force-strain curve by Thelen. Obtained in Matlab 26 Figure 3.9. Tendon force-strain curve by Thelen. Obtained in Matlab 27 Figure 3.10. Active force-length curve by Thelen. Obtained in Matlab 27 Figure 3.11. Force velocity curve by Thelen. Obtained in Matlab 28 Figure 3.12. Tendon force-length curve by De Groote-Fregly. Obtained in Matlab 29 Figure 3.13. Passive force-length curve by De Groote-Fregly. Obtained in Matlab 29 Figure 3.6. Active force length curve. OpenSim Documentation [28] 24 Figure 3.7. Tendon force length curve. OpenSim Documentation [28] 25 Figure 3.8. Tendon force-strain curve by Thelen. Obtained in Matlab 26 Figure 3.9. Tendon force-strain curve by Thelen. Obtained in Matlab 27 Figure 3.10. Active force-length curve by Thelen. Obtained in Matlab 27 Figure 3.11. Force velocity curve by Thelen. Obtained in Matlab 28 Figure 3.12. Tendon force-length curve by De Groote-Fregly. Obtained in Matlab 29 Figure 3.13. Passive force-length curve by De Groote-Fregly. Obtained in Matlab 30 Figure 3.14. Active force-length curve by De Groote-Fregly. Obtained in Matlab 31 Figure 3.15. Force-velocity curve by De Groote-Fregly. Obtained in Matlab. 32 Figure 4.1. Example of open-loop mechanism in OpenSim GUI: 3-bar pendulum. [6] 36 Figure 4.2. Operating scheme of Inverse Kinematics in OpenSim 37 Figure 4.3. Operating scheme of Inverse Dynamics in OpenSim 38 Figure 4.4. Example of closed-loop model in OpenSim [31] 38 Figure 5.1. 4-bar linkage in OpenSim. [6] 43 Figure 5.2. The 4-bar at t = 0 46 Figure 5.3. The 4-bar at t = 1.5 s 46 Figure 5.4. The 4-bar at t = 3 s 46 Figure 5.5. The 4-bar at t = 4.5 s 46 Figure 5.6. The 4-bar at t = 6 s 46 Figure 5.7. 4-bar mechanism built in WinMecC with generalized coordinates indicated [34] 47 Figure 5.8. 4-bar linkage with two muscles, each at one side. OpenSim [6] 53 Figure 6.1. Analytical motor torque of 4-bar. Obtained in Matlab 58 Figure 6.2. Analytical power of 4-bar. Obtained in Matlab 58 Figure 6.3. Power resulting in ID by OpenSim. Obtained in Matlab 59 Figure 6.4. Analytical power of 4-bar. Obtained in Matlab 60 Figure 6.5. Power in a non-muscled 4-bar mechanism. Obtained in Matlab 60 Figure 6.6. Comparison of total power obtained in Moco and the analytical way. Obtained in Matlab 60 Figure 6.7. Power resulted in muscled 4-bar. Obtained in Matlab. 61 Figure 6.8. Comparison of total power obtained in Moco and the analytical way, in the muscled 4-bar. Obtained in Matlab 62 Figure 6.9. Contribution of muscles and Coord. Actuators to power. Obtained in Matlab 62 Figure 6.10. Comparison of powers by muscle and their activations. Obtained in Matlab 63 Figure 6.11. Right muscle forces. Obtained in Matlab 64 Figure 6.12. Left muscle forces. Obtained in Matlab 64 Figures 6.13. New insertion point of muscles in input bar. OpenSim [6] 65 Figures 6.14. New insertion point of muscles in input bar. OpenSim [6] 65 Figure 6.15. 4-muscles 4-bar linkage in OpenSim [6] 65 Figure 6.16. Influence of insertion point in muscle power. Obtained by Matlab 66 Figure 6.17. Influence of insertion point in actuators’ power. Obtained by Matlab 67 Figure 9.1. Normalized tendon force in both muscles. OpenSim Moco [13] 76 Figure 9.2. Normalized length force in both muscles. OpenSim Moco [13] 76 Figure 9.3. Normalized tendon force in right muscle. OpenSim Moco [13] 77 Figure 9.4. Normalized tendon force in left muscle. OpenSim Moco [13] 77 Figure 9.5. Muscle lengths. OpenSim Moco [13] 78 Figure 9.6. Muscle fiber lengths. OpenSim Moco [13] 78 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons xxi Figure 9.7. Normalized contraction velocities in both muscles. OpenSim Moco [13] 79 Figure 9.8. Normalized passive force in both models, right muscle. OpenSim Moco [13] 79 Figure 9.9. Normalized left passive force in both models. OpenSim Moco [13] 80 Figure 9.10. Right active force-length multiplier in both models. OpenSim Moco [13] 80 Figure 9.11. Left active force-length multiplier in both models. OpenSim Moco [13] 81 Figure 9.12. Right active force-velocity multiplier in both models. OpenSim Moco [13] 81 Figure 9.13. Left active force-velocity multiplier in both models. OpenSim Moco [13] 82 Figure 9.14. Normalized active forces in both muscles. OpenSim Moco [13] 82 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons xxiii 1 1 INTRODUCTION The implementation of exoskeletons in human people is one of the most complex and exciting challenges being faced by Biomechanics since the last few decades. These devices, whose origin can be found on the military field [2], help people to perform diverse types of motion by supporting an extra amount of energy. There are two types of exoskeletons according to the source of energy [1]: active exoskeletons, if they contain an electrical device which supports energy from a battery or other sort of external source like pneumatical or hydraulic circuits; and passive exoskeletons, whose way to obtain energy consists of the springs, dampers and other elastic elements which store energy and release it during the motion. The biomechanical models used in this work simulate passive exoskeletons, because no electrical actuator is implemented. Figures 1.1 and 1.2: Industrial worker wearing an Exhauss passive exoskeleton (Theurel et al., [3]) and an active exoskeleton worn in lower limbs (Yale University, [4]) The study of dynamics in the whole human-exoskeleton is not an easy problem to solve, as two complex systems are being considered, in the sense of the great magnitude of bodies, joints and muscles they present. Moreover, the only form to attach an exoskeleton is by linking it to the person through two points, creating therefore a closed loop which complicates even more the problem. This project has the goal to evaluate the adequacy of OpenSim Moco (Stanford University, [5]) to solve Inverse Dynamics problems (it is said, problems where forces and moments are calculated parting from a given and perfectly defined kinematics) in musculoskeletal models with exoskeletons by minimizing a determined cost function, it is said, solving an optimal control problem where the equations of motion are part of the constraints in the system. OpenSim Moco is a software developed by Stanford University [5] which solves optimal control problems by direct collocation methods that will be detailed in Chapter 3. Moco is part of OpenSim, an opensource software system also powered by Stanford University [6], which serves to create musculoskeletal models and perform simulations with them. Every model in this work have been built with OpenSim. The difficult to study such a complex closed-loop model in OpenSim like a human wearing an exoskeleton has OpenSim Moco 8 8 Now it is time to describe the methods used by Moco to solve the optimal control problem: direct collocation methods. 2.3 Direct collocation methods The problem just presented, the optimal control, is a trajectory optimization problem continuous in time, as it introduces integrals in the objective function. To solve it, Moco makes use of direct collocation methods [13]. These algorithms transform the optimal control problem into a nonlinear discrete optimization problem, by discretizing the trajectory optimization problem into a determined number of control points, named nodes [15]. The continuous functions belonging to the original optimal control problem are substituted by polynomial splines whose value in control points must be the same as in the continuous functions. A graphic representation is shown in Figure 2.2: Figure 2.2. Lineal interpolation through nodes in optimal control problem. Kelly, M. [15] One of the main advantages about using direct collocation resides on the fact that the dynamic constraints (it is said, the equations of motion), are taken in count in the cost function, as they include control variables which must minimize the goal [12, 15]. This makes slightly easier to find a solution for the optimization problem. Moco can use two different collocation methods to solve the optimal control problem: the trapezoidal method, and the Hermite-Simpson’s method [13]. 2.3.1 Trapezoidal collocation The trapezoidal direct collocation method uses the trapezoidal quadrature rule to discretize the continuous problem in N control points and approximate the integrals as summations of discrete variables, evaluated in those control nodes [15]: ∫𝑓(𝜏,𝑥(𝜏)) 𝑑𝜏 𝑡𝑓 𝑡0≈ ∑𝑡𝑟𝑎𝑝𝑖=∑1 2∗(𝑡𝑖+1−𝑡𝑖)∗(𝑓𝑖+𝑓𝑖+1) 𝑁−1 𝑖=0 𝑁−1 𝑖=0 (2.20) Where N is the total number of control points and 𝑓(𝜏,𝑥(𝜏)) is any function belonging to an integrand in the optimization problem. It is needed that the following condition gets respected in each control point [15]: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 9 𝑓(𝑡𝑖)=𝑓𝑖  (2.21) It is said, in control points, the value of the continuous function and its approximative spline must be the same. The process to obtain the values of control and state variables can be explained as follows in the next lines. In first place, the spline used to approximate the control trajectories is a linear polynomial [15]: 𝑢(𝑡)≈𝑢𝑘+𝜏 ℎ𝑘∗(𝑢𝑘+𝑢𝑘+1) (2.22) Where 𝜏 and ℎ𝑘 are expressed as follows: 𝜏=𝑡−𝑡𝑘 (2.23) ℎ𝑘=𝑡𝑘+1−𝑡𝑘 (2.24) Regarding states, the dynamic equations can be represented in function of the first derivative of the state variables: 𝑥󰇗=𝑓 (2.25) If trapezoidal quadrature rule is used, the following relation is obtained [15]: 𝑓≈𝑓𝑘+𝜏 ℎ𝑘∗(𝑓𝑘+𝑓𝑘+1) (2.26) State trajectories are obtained by integrating this equation: 𝑥≈𝑥𝑘+𝑓𝑘∗𝑡+𝜏2 2ℎ𝑘∗(𝑓𝑘+𝑓𝑘+1) (2.27) This way, after applying the trapezoidal rule, the following non-linear discrete optimization problem appears [13]: min∑𝑤𝑗𝐽𝑗(𝑡0,𝑡𝑓, 𝑗𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑐,𝑗)+𝑤𝜆∗∑𝑡𝑟𝑎𝑝𝑖( 𝑛 𝑖=1 ‖𝜆𝑖‖2 2) 𝑆𝑐,𝑗=∑𝑡𝑟𝑎𝑝𝑖(𝑠𝑐,𝑗(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑁 𝑖=0 OpenSim Moco 10 10 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑢𝑖=𝑢𝑖−1+𝑡𝑟𝑎𝑝𝑖(𝑣) 𝑖=1,…,𝑁 𝑀(𝑞𝑖,𝑝)∗𝑣𝑖+𝐺(𝑞,𝑝)𝑇∗𝜆𝑖=𝑓𝑎𝑝𝑝(𝑡𝑖,𝑦𝑖,𝑥𝑖,𝑝)−𝑓𝑖𝑛𝑒𝑟𝑡𝑖𝑎𝑙(𝑞𝑖,𝑢𝑖,𝑝) 𝑖=0,…,𝑁 𝑧𝑒𝑥𝑖=𝑧𝑒𝑥,𝑖−1+𝑡𝑟𝑎𝑝𝑖(𝑓𝑒𝑥,𝑧󰇗(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑖=1,…,𝑁 𝑧𝑖𝑚𝑖=𝑧𝑖𝑚,𝑖−1+𝑡𝑟𝑎𝑝𝑖(𝜁) 𝑖=1,…,𝑁 0=𝑓𝑧󰇗𝑖𝑚(𝑡)(𝑡𝑖,𝑦𝑖,𝜁𝑖,𝑥𝑖,𝜆𝑖,𝑝) 𝑖=0,…,𝑁 𝛷(𝑞𝑖,𝑝)=0 𝑖=0,…,𝑁 𝑉𝐿,𝑘≤𝑉𝑘(𝑡0,𝑡𝑓,𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑏,𝑘)≤𝑉𝑈,𝑘 𝑆𝑏,𝑘=∑𝑡𝑟𝑎𝑝𝑖(𝑠𝑏,𝑘(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑁 𝑖=0 𝑘=1,…,𝐾 𝑔𝐿≤𝑔(𝑡𝑖,𝑦𝑖,𝑥𝑖,𝜆𝑖,𝑝)≤𝑔𝑈 𝑖=1,…,𝑁 𝑤𝑖𝑡ℎ 𝑟𝑒𝑠𝑝𝑒𝑐𝑡 𝑡𝑜 𝑣𝑖∈[−𝑣𝐵,𝑣𝐵] 𝑦0,𝐿≤𝑦0≤𝑦0,𝑈 𝑦𝑓,𝐿≤𝑦𝑛≤𝑦𝑓,𝑈 𝑦𝐿≤𝑦𝑖≤𝑦𝑈 𝑖=0,…,𝑁−1 𝑥0,𝐿≤𝑥0≤𝑥0,𝑈 𝑥𝑓,𝐿≤𝑥𝑛≤𝑥𝑓,𝑈 𝑥𝐿≤𝑥𝑖≤𝑥𝑈 𝑖=0,…,𝑁−1 𝑡0,𝐿≤𝑡0≤𝑡0,𝑈 𝑡𝑓,𝐿≤𝑡𝑛≤𝑡𝑓,𝑈 𝑦𝐿≤𝑦(𝑡)≤𝑦𝑈 𝑥𝐿≤𝑥(𝑡)≤𝑥𝑈 𝑝𝐿≤𝑝≤𝑝𝑈 𝜁𝐿≤𝜁𝑖≤𝜁𝑈 𝑖=0,…,𝑁 𝜆𝐿≤𝜆𝑖≤𝜆𝑈 𝑖=0,…,𝑁 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 11 Where 𝑣 is generalized acceleration and 𝑣𝐵 a very large quantity [13]. 2.3.2 Hermite-Simpson collocation This algorithm is the one used by default by Moco. It allows more accurate results than the trapezoidal rule because it introduces additional collocation points at the mesh interval midpoints, resulting to a total of 2N+1 grid points in which the trajectory is discretized [15]. This way, control trajectories are approximated by quadratic splines, while state trajectories are approximated to cubic splines. First of all, it is necessary to define the Simpson’s rule of quadrature [15]: ∫𝑓(𝜏,𝑥(𝜏)) 𝑑𝜏 𝑡𝑓 𝑡0≈ ∑𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖=∑ℎ𝑖 6∗(𝑓𝑖−1+4∗𝑓𝑖+𝑓𝑖+1) 𝑁−1 𝑖=0 𝑁−1 𝑖=0 (2.28) With ℎ𝑖 defined as in Equation. The value of the state variable in the midpoint belonging to the instant tk+1/2 is obtained with the Hermite interpolant: ℎ𝑒𝑟𝑚𝑖𝑡𝑒𝑘+1/2=𝑥𝑘+1 2=1 2∗(𝑥𝑘+𝑥𝑘+1)+ℎ𝑘 8∗(𝑓𝑘−𝑓𝑘+1) (2.29) Following equations to obtain the control and state variables through time are demonstrated in Kelly, M. [15]: 𝑢(𝑡)=2 ℎ𝑘 2∗(𝜏−ℎ𝑘 2)∗(𝜏−ℎ𝑘)∗𝑢𝑘−4 ℎ𝑘 2∗(𝜏)∗(𝜏−ℎ𝑘)∗𝑢𝑘+1 2+2 ℎ𝑘 2∗(𝜏)∗(𝜏−ℎ𝑘)∗𝑢𝑘+1 (2.30) 𝑥(𝑡)=𝑥𝑘+𝑓𝑘(𝜏 ℎ𝑘)+1 2∗(−3𝑓𝑘+4𝑓𝑘+1 2−𝑓𝑘+1)∗(𝜏 ℎ𝑘)2+1 3∗(2𝑓𝑘−4𝑓𝑘+1 2+2𝑓𝑘+1)∗(𝜏 ℎ𝑘)3(2.31) Therefore, the formulation of the optimal control problem after applying the Hermite-Simpson collocation method would be the following one [13]: min∑𝑤𝑗𝐽𝑗(𝑡0,𝑡𝑓, 𝑗𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑐,𝑗)+𝑤𝜆∗∑𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖( 𝑛 𝑖=1 ‖𝜆𝑖‖2 2) 𝑆𝑐,𝑗=∑𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖(𝑠𝑐,𝑗(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑁 𝑖=0 OpenSim Moco 12 12 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑢𝑖=ℎ𝑒𝑟𝑚𝑖𝑡𝑒𝑖(𝑢,𝑣) 𝑖=1…𝑁 𝑢𝑖=𝑢𝑖−1+𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖(𝑣) 𝑖=1…𝑁 𝑀(𝑞𝑖,𝑝)∗𝑣𝑖+𝐺(𝑞𝑖,𝑝)𝑇∗𝜆𝑖=𝑓𝑎𝑝𝑝(𝑡𝑖,𝑦𝑖,𝑥𝑖,𝑝)−𝑓𝑖𝑛𝑒𝑟𝑡𝑖𝑎𝑙(𝑞𝑖,𝑢𝑖,𝑝) 𝑖=0,…,𝑁 𝑀(𝑞𝑖 ,𝑝)∗𝑣𝑖+𝐺(𝑞𝑖 ,𝑝)𝑇∗𝜆𝑖=𝑓𝑎𝑝𝑝(𝑡𝑖 ,𝑦𝑖 ,𝑥𝑖,𝑝)−𝑓𝑖𝑛𝑒𝑟𝑡𝑖𝑎𝑙(𝑞𝑖 ,𝑢𝑖 ,𝑝) 𝑖=1,…,𝑁 𝑧𝑒𝑥𝑖=𝑧𝑒𝑥,𝑖−1+𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖(𝑓𝑒𝑥,𝑧󰇗(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑖=1,…,𝑁 𝑧𝑒𝑥𝑖=ℎ𝑒𝑟𝑚𝑖𝑡𝑒𝑖(𝑓𝑒𝑥,𝑧󰇗(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑖=1,…,𝑁 𝑧𝑖𝑚𝑖=ℎ𝑒𝑟𝑚𝑖𝑡𝑒𝑖(𝜁) 𝑖=1,…,𝑁 𝑧𝑖𝑚𝑖=𝑧𝑖𝑚,𝑖−1+𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖(𝜁) 𝑖=1,…,𝑁 0=𝑓𝑧󰇗𝑖𝑚(𝑡)(𝑡𝑖,𝑦𝑖,𝜁𝑖,𝑥𝑖,𝜆𝑖,𝑝) 𝑖=0,…,𝑁 𝑥𝑖 =𝑥𝑖−1+𝑥𝑖 2 𝑖=1,…,𝑁 𝛷(𝑞𝑖,𝑝)=0 𝑖=0,…,𝑁 𝛷󰇗(𝑞𝑖,𝑢𝑖,𝑝)=0 𝑖=0,…,𝑁 𝛷󰇘(𝑡𝑖,𝑦𝑖,𝑥𝑖,𝜆𝑖,𝑝)=0 𝑖=0,…,𝑁 𝑉𝐿,𝑘≤𝑉𝑘(𝑡0,𝑡𝑓,𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑏,𝑘)≤𝑉𝑈,𝑘 𝑆𝑏,𝑘=∑𝑠𝑖𝑚𝑝𝑠𝑜𝑛𝑖(𝑠𝑏,𝑘(𝑡,𝑦,𝑥,𝜆,𝑝)) 𝑘=1,…,𝐾 𝑁 𝑖=0 𝑔𝐿≤𝑔(𝑡𝑖,𝑦𝑖,𝑥𝑖,𝜆𝑖,𝑝)≤𝑔𝑈 𝑖=1,…,𝑁 𝑤𝑖𝑡ℎ 𝑟𝑒𝑠𝑝𝑒𝑐𝑡 𝑡𝑜 𝑦0,𝐿≤𝑦0≤𝑦0,𝑈 𝑦𝑓,𝐿≤𝑦𝑛≤𝑦𝑓,𝑈 𝑦𝐿≤𝑦𝑖≤𝑦𝑈 𝑖=0,…,𝑁−1 𝑦𝐿≤𝑦𝑖≤𝑦𝑈 𝑖=0,…,𝑁 𝑥0,𝐿≤𝑥0≤𝑥0,𝑈 𝑥𝑓,𝐿≤𝑥𝑛≤𝑥𝑓,𝑈 𝑥𝐿≤𝑥𝑖≤𝑥𝑈 𝑖=0,…,𝑁−1 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 13 𝑥𝐿≤𝑥𝑖≤𝑥𝑈 𝑖=0,…,𝑁 𝑡0,𝐿≤𝑡0≤𝑡0,𝑈 𝑡𝑓,𝐿≤𝑡𝑛≤𝑡𝑓,𝑈 𝑝𝐿≤𝑝≤𝑝𝑈 𝜁𝐿≤𝜁𝑖≤𝜁𝑈 𝑖=0,…,𝑁 𝜁𝐿≤𝜁𝑖≤𝜁𝑈 𝑖=0,…,𝑁 𝜆𝐿≤𝜆𝑖≤𝜆𝑈 𝑖=0,…,𝑁 𝜆𝐿≤𝜆𝑖≤𝜆𝑈 𝑖=0,…,𝑁 2.4 Types of Moco problems Moco can solve fundamentally three types of problems, depending on the presence or absence of information about the model’s motion [12, 13]. 2.4.1 MocoInverse This is the type of problem faced in this work. In MocoInverse, model’s kinematics is completely prescribed, so the objective is to obtain the forces which produce this motion, via the control variables present in the objective function, and the actuator states [12, 13]. Therefore, it deals with an Inverse Dynamic problem. The scheme in Figure helps to understand it: Figure 2.3. Operating scheme of MocoInverse problems This is the type of problem approached in this work, so it is useful to know the formulation of the problem. It is presented right below: MOCO OpenSim model Cost function Constraints States Controls Motion OpenSim Moco 14 14 min∑𝑤𝑗𝐽𝑗(𝑡0,𝑡𝑓, 𝑗𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑐,𝑗) 𝑆𝑐,𝑗=∫ 𝑠𝑐,𝑗(𝑡,𝑦,𝑥,𝜆,𝑝)𝑑𝑡 𝑡𝑓 𝑡0 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑢=𝑞󰇗 𝑀(𝑞,𝑝)∗𝑢󰇗+𝐺(𝑞,𝑝)𝑇∗𝜆=𝑓𝑎𝑝𝑝(𝑡,𝑦,𝑥,𝑝)−𝑓𝑖𝑛𝑒𝑟𝑡𝑖𝑎𝑙(𝑞,𝑢,𝑝) 𝑧󰇗𝑒𝑥(𝑡)=𝑓𝑧󰇗𝑒𝑥(𝑡)(𝑡,𝑦,𝑥,𝜆,𝑝) 0=𝑓𝑧󰇗𝑖𝑚(𝑡)(𝑡,𝑦,𝑧󰇗𝑖𝑚(𝑡),𝑥,𝜆,𝑝) 𝑉𝐿,𝑘≤𝑉𝑘(𝑡0,𝑡𝑓,𝑦0,𝑦𝑓,𝜆0,𝜆𝑓,𝑥0,𝑥𝑓,𝑝,𝑆𝑏,𝑘)≤𝑉𝑈,𝑘 𝑆𝑏,𝑘=∫ 𝑠𝑏,𝑘(𝑡,𝑦,𝑥,𝜆,𝑝)𝑑𝑡 𝑡𝑓 𝑡0 𝑘=1,…,𝐾 𝑔𝐿≤𝑔(𝑡,𝑦,𝑥,𝜆,𝑝)≤𝑔𝑈 𝑤𝑖𝑡ℎ 𝑟𝑒𝑠𝑝𝑒𝑐𝑡 𝑡𝑜 𝑦0,𝐿≤𝑦0≤𝑦0,𝑈 𝑦𝑓,𝐿≤𝑦𝑓≤𝑦𝑓,𝑈 𝑥0,𝐿≤𝑥0≤𝑥0,𝑈 𝑥𝑓,𝐿≤𝑥𝑓≤𝑥𝑓,𝑈 𝑡0,𝐿≤𝑡0≤𝑡0,𝑈 𝑡𝑓,𝐿≤𝑡𝑓≤𝑡𝑓,𝑈 𝑦𝐿≤𝑦(𝑡)≤𝑦𝑈 𝑥𝐿≤𝑥(𝑡)≤𝑥𝑈 𝑝𝐿≤𝑝≤𝑝𝑈 𝜆(𝑡) Computation time in MocoInverse is shorter than in the other problems due to the low number of unknowns, therefore it is the simplest problem which can be faced with Moco [12, 13]. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 15 2.4.2 MocoTrack In these types of problem, there exists information about the motion performed by the model, this information normally consists of experimental data captured by sensors. However, kinematics becomes unknown, so it is one of the variables to obtain, along with controls and states. The objective deals with fulfill the equations of motion (and the rest of constraints) while minimizing as much as possible the deviation of the resulting motion from the experimental data [12, 13]. Figure shows the representative scheme of problems of this nature: Figure 2.4. Operating scheme in MocoTrack problems It takes longer to compute in comparison to MocoInverse, as the motion is part of the objective function [12]. 2.4.3 MocoStudy It is the most complex case study approached in Moco, and the most general at the same time, as its formulation is equivalent to the one presented in the sub-chapter dedicated to the optimal control problem. Motion is completely unknown, as there is no data about it, therefore, its computation time is the longest by far [12, 13]. This is the scheme corresponding to MocoStudy: MOCO OpenSim model Motion data (experimental) Motion Controls States Cost function Constraints OpenSim Moco 16 16 Figure 2.5. Operating scheme of MocoStudy problems. MOCO OpenSim model Motion Controls Parameters Cost function Constraints Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 17 Muscle-tendon modelling 24 24 ● Active force length curve: the optimal fiber length just described is that one which allows the muscle to generate the maximum force possible. If it is surpassed, not only passive forces start to appear, but also active forces get reduced. Of course, under this threshold active forces decrease too. Diverse representations of the curve may be appreciated in Figure 3.6: Figure 3.6. Active force length curve. OpenSim Documentation [28] ● Force velocity curve: in the movement, a muscle can be in two different phases: lengthening and shortening. In the first phase, which corresponds to a negative contraction velocity (v < 0) the muscle decreases the generated force, while the second phase deals with positive contraction velocity (v > 0). The muscles can generate more force than in the previous phase. The Figure 3.7 shows some representations of these curve. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 25 Figure 3.7. Tendon force length curve. OpenSim Documentation [28] Now that the process to obtain different forces has been explained, some muscle models will be presented. 3.4 Muscle models 3.4.1 Thelen muscle model This model was developed by Thelen [18]. Starting with the tendon force, in Thelen’s model it is necessary to define, at first, tendon strain, ɛt, as the relative displacement of the tendon with respect to its slack length: 𝜀𝑡=𝑙𝑠𝑒−𝑙𝑠𝑙𝑎𝑐𝑘 𝑙𝑠𝑙𝑎𝑐𝑘 (3.20) Now the force exerted by the tendon according to Thelen can be expressed: 𝐹𝑡=𝐹𝑚𝑎𝑥∗0.10377∗(𝑒91∗𝜀𝑡−1) 𝑖𝑓 0<𝜀𝑡<0.01516 (3.21) 𝐹𝑡=𝐹𝑚𝑎𝑥∗(37.526∗𝜀𝑡−0.26029) 𝑖𝑓 𝜀𝑡 ≥0.01516 (3.22) The Figure 3.8 shows the graphic representation of this curve. The exponential stretch is appreciated from the beginning until the strain reaches a value of 𝜀𝑡= 0.01516, where the lineal part starts. Figure 3.8. Tendon force-strain curve by Thelen. Obtained in Matlab. With respect to the parallel element, which produces passive forces, the expression depends on the normalized Muscle-tendon modelling 26 26 fiber length, obtained as follows: 𝑙𝑝𝑒 =𝑙𝑝𝑒 𝑙𝑜𝑝𝑡 (3.23) Being 𝑙𝑜𝑝𝑡 the optimal fiber length. This way, the equation to obtain the normalized passive force is the following one: 𝑓𝑝𝑒 =1+𝑘𝑝𝑒 𝜀𝑐𝑒 ∗(𝑙󰆻𝑐𝑒−(1+𝜀𝑐𝑒)) 𝑖𝑓 𝑙󰆻𝑐𝑒> 1+𝜀𝑐𝑒 (3.24) 𝑓𝑝𝑒 =𝑒𝑘𝑝∗(𝑙𝑐𝑒 −1) 𝜀𝑐𝑒 𝑒𝑘𝑝𝑒 𝑖𝑓 𝑙󰆻𝑐𝑒≤ 1+𝜀𝑐𝑒 (3.25) Where 𝑘𝑝𝑒 is a shape factor got experimentally, and 𝜀𝑐𝑒 is the passive muscle strain at the maximum isometrical force (in this work, 𝑘𝑝𝑒 = 4 and 𝜀𝑐𝑒 = 0.6). The graphic representation comes as follows in the Figure 3.9: Figure 3.9. Tendon force-strain curve by Thelen. Obtained in Matlab. The next step would be to calculate the effect of the fiber length in the active force. So that, the next exponential expression is used: 𝑓𝐿=𝑒−(𝑙󰆻𝑐𝑒−1)2 𝛾 (3.26) Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 27 Where 𝛾 corresponds to the half width of the curve 1 𝑒 . In this project, a value of 0.45 has been assigned to this parameter. The Figure 3.10 shows the graphic relation between active force length parameter and the normalized fiber length: Figure 3.10. Active force-length curve by Thelen. Obtained in Matlab Last, the relationship force-velocity must be represented. This expression is written in function of the normalized contraction velocity, which comes defined as follows: 𝑣𝑐𝑒 = 𝑣𝑐𝑒 𝑣𝑚𝑎𝑥 (3.27) Resulting 𝑣𝑚𝑎𝑥 the maximum contraction velocity of the muscle, in this work it is assigned a 10 m/s value. Now the expression force-velocity can be written: 𝑓𝑣=1+𝑣𝑐𝑒 1−𝑣 𝑐𝑒 𝐴𝑓 𝑖𝑓 𝑣𝑐𝑒 ≤ 0 (3.28) 𝑓𝑣=1+𝑣𝑐𝑒∗𝑓𝑣𝑚𝑎𝑥∗(2+2 𝐴𝑓) 1+𝑣𝑐𝑒∗(2+2 𝐴𝑓) 𝑓𝑣𝑚𝑎𝑥−1 𝑖𝑓 𝑣𝑐𝑒>0 (3.29) Where 𝑓𝑣𝑚𝑎𝑥 is the maximum normalized lengthening force (its assigned value is 1.8), and 𝐴𝑓 a shape factor (its value has been set to 0.25). The curve representing these equations is displayed in Figure 3.11: Muscle-tendon modelling 28 28 Figure 3.11. Force velocity curve by Thelen. Obtained in Matlab 3.4.2 De Groote-Fregly muscle model This model was developed by De Groote, Fregly et al. [19]. Unlike Thelen’s model, which defined a tendon strain, De Groote-Fregly’s considers a normalized tendon length as follows: 𝑙󰆻𝑠𝑒=𝑙𝑠𝑒 𝑙𝑠𝑙𝑎𝑐𝑘 (3.29) Being 𝑙𝑠𝑙𝑎𝑐𝑘 the tendon slack length. This adimensional magnitude has a lower bound of zero, as the tendon cannot have a shorter length than its slack length. Known this, the equation relating the tendon force and its length is written the following way: 𝑙󰆻𝑠𝑒=𝑙𝑛 (5∗𝑓󰆻𝑠𝑒+0.25) 48 +0.995 (3.30) Its graphic representation comes displayed as follows in Figure 3.12: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 29 Figure 3.12. Tendon force-length curve by De Groote-Fregly. Obtained in Matlab Concerning the passive force-length relationship, De Groote-Fregly’s model utilizes the same expression as Thelen’s, but without considering the lineal stretch, the whole curve corresponds to the exponential expression: 𝑓𝑝𝑒 =𝑒𝑘𝑝∗(𝑙𝑐𝑒 −1) 𝜀𝑐𝑒 𝑒𝑘𝑝𝑒 (3.31) Where, as already explained in Thelen’s model, 𝑘𝑝 corresponds to an experimental shape factor (a value of 4 has been stablished for it), while 𝜀𝑐𝑒 is the passive muscle strain at maximum isometric muscle force (it has been set to 0.6). The curve corresponding to the equation is shown in Figure 3.13: Figure 3.13. Passive force-length curve by De Groote-Fregly. Obtained in Matlab. Muscle-tendon modelling 30 30 With regards to the active force length characteristic curve, De Groote-Fregly’s model sets it as the sum of three Gaussian expressions, which contain several experimental values. These relationships are written as follows: 𝑓𝑙= ∑𝑏1𝑖∗𝑒(−0.5∗(𝑙𝑐𝑒 −𝑏2𝑖)2 𝑏3𝑖+𝑏4𝑖∗𝑙𝑐𝑒 ) 3 𝑖=1 (3.32) Being bij, i=1…4, j=1…3 experimental parameters whose value comes collected in the Table 3.1: Table 3.1. Experimental parameters in passive force-length equation by De Groote-Fregly [19] b11 b21 b31 b41 b12 b22 b32 b42 b13 b23 b33 b43 0.815 1.055 0.162 0.063 0.433 0.717 -0.03 0.02 0.1 1 0.354 0 Its graphic representation comes displayed in the Figure 3.14: Figure 3.14. Active force-length curve by De Groote-Fregly. Obtained in Matlab With respect to the force velocity curve, it is described by a logarithmic function which contains some predefined parameters: 𝑓𝑣=𝑑1∗𝑙𝑛 [(𝑑2∗𝑣𝑐𝑒 +𝑑3)+√((𝑑2∗𝑣𝑐𝑒 +𝑑3)2+1)] +𝑑4 (3.33) Where the parameters di , i = 1…4 have been assigned the next values: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 31 Table 3.2. Parameters in force-velocity equation by De Groote-Fregly [19] 𝑑1 𝑑2 𝑑3 𝑑4 -0.3183 -8.1492 -0.3741 0.8856 The force-velocity curve according to De Groote-Fregly’s model results to be the following one: Figure 3.15. Force-velocity curve by De Groote-Fregly. Obtained in Matlab. Muscle-tendon modelling 32 32 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 33 4 THE CLOSED-LOOP PROBLEM As already explained, the intention of this work is to study the possibility of using OpenSim Moco to study dynamics in human people wearing exoskeletons. When attached to a person, the union formed by the individual and the exoskeleton constitutes a muscled closed-loop mechanism. This way, the objective of the project is summed up in checking if OpenSim Moco can solve Inverse Dynamics problems, firstly, in closed-loop mechanisms, and then, in muscled closed-loop mechanisms. This chapter explains, in a theoretical and concise manner, how the objective right just described is faced in this work. It starts with the way OpenSim approaches Inverse Dynamics, as well as its form to solve open-loop and closed-loop mechanisms, and the use of Moco to face these problems. The four-bar linkage, the closed-chain mechanism used in this project, will be presented. By last, the analytical inverse dynamics problem applied to a multibody system will be introduced, as this is the problem whose result will be compared to the Inverse Dynamics made in Moco. 4.1 Open loop mechanisms in OpenSim As introduced earlier, OpenSim allows to develop musculoskeletal models and dynamic simulations to study its movement. Equations of motion are formulated by OpenSim by using Simbody, an open-source multibody dynamics solver [20]. OpenSim models are constituted by bodies, joints and, if needed, constraints: a) Bodies: they are rigid solids connected by joints, and conform the main blocks in the model. Generally, each body has assigned a joint which links it to a previous body (this body is considered a parent one), and it is also connected to a successive body (a child one) by another joint [20]. b) Joints: they serve to stablish the unions between bodies. Models are built following a chain of bodies linked by joints, where each joint relates two bodies: the parent and the child one, except for the initial body (normally the ground, it is always a parent one) and the final bodies (there are as many as different chains the model contains; final bodies are always child ones).[20] Joints have a great importance as they define the coordinates of the model, it is said, the parameters which will determine the motion of the model. In building the joint, not only the coordinates must be declared, but also the position and orientation with respect to the parent and child bodies. There are different available joints in OpenSim according to the number of coordinates they introduce, here the most used in this Project will be described: ● Ball joints: they introduce 3 rotational coordinates between bodies, blocking any possible translation between them. ● Pin joints: they only introduce one rotational coordinate, through the Z-axis common to both bodies, restricting the other two restrictions and the three translations. ● Weld joints: they fuse the two bodies, avoiding any possible movement, rotational or translational. No coordinates are defined in Weld joints. The Closed-loop problem 40 40 As the accelerations are already known, Lagrange multipliers can be easily achieved: 𝝀=𝑪𝒒 𝑻−𝟏(𝑸𝒂𝒑𝒍+𝑸𝒗−𝑴∗𝒒󰇘) (4.9) Now, generalized forces can be calculated: 𝑸𝒈𝒆𝒏=−𝑪𝒒 𝑻∗𝝀 (4.10) Ending with the Inverse Dynamics problem. These generalized forces represent the forces required to apply in each generalized coordinate to fulfill each restriction; in fact, there is a generalized coordinate associated to each generalized force, and one generalized force associated to each restriction. In the next chapter, the whole theoretical set just presented will be particularized for the 4-bar linkage. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 41 Inverse Dynamics in the 4-bar linkage 42 42 5 INVERSE DYNAMICS IN THE 4-BAR LINKAGE In this chapter, the 4-bar mechanism object of study in this project will be presented, including the description of each bar and the generalized coordinates in the problem, as well as the subsequent modifications made later (the inclusion of two muscles to study the interaction of the muscles with the bars of the mechanism) and the motion subject of study the mechanism is doing. Then, the Inverse Dynamics problem particularized for the 4-bar mechanism will be introduced and slightly detailed, for both the analytical and the MocoInverse problem, presenting all variables, equations and formulations needed to solve it. In Chapter 7, results from all these problems will be displayed. 5.1 The 4-bar linkage 5.1.1 Introduction and composition The four-bar linkage is the simplest closed chain mechanism which can be built. This mechanism has been used to prove the validity of the results in a closed chain mechanism, in this case, the simplest possible, as it only has one degree of freedom. Figure 5.1. 4-bar linkage in OpenSim. [6] Figure 5.1 shows the four bar linkage used in this project. Different bars composing the mechanism are the following ones: a) Fixed bar: it does not move as it matches the ground. In this case, it is 1.2 meters long. b) Input bar: it is linked to the fixed bar through a ball joint (later the reason not to be a pin joint will be explained), and it performs a 360 degrees spin through the full motion cycle. Its length is 0.4 m. c) Coupler bar: it is connected to the input bar through a pin joint, and transmits the motion to the output bar. It is 1 m long. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 43 d) Output bar: it receives the motion from the coupler bar, whereby it is connected through a pin joint. It is 0.8 m long. In the literature [32], the output bar is connected to the ground by a pin joint, but as previously explained, in OpenSim a chain cannot be closed only by implementing joints, constraints must be used on their behalf. The proceed followed to close the mechanism here consisted of create a new bar (the one laid in the ground in Figure), and link it to the output bar through a pin joint, and to the fixed bar by a weld constraint. This is equivalent to join the mechanism to the ground with a pin joint. 5.1.2 Grashof condition In order to perform a complete 360 degrees revolution, the four-bar mechanism must fulfill the Grashof condition, which states that the sum of the lengths of the shortest and the longest bar must be less or equal than the sum of the lengths of the other two bars [32]. This come expressed in the Equation 𝐿1+𝐿4≤𝐿2+𝐿3 (5.1) Being L1 and L4 the longest and shortest bars respectively, and L2 and L3 the other two. In this case, L1 = 1.2 m, L2 = 1 m, L3 = 0.8 m and L4 = 0.4 m, so the Grashof condition gets fulfilled: 1.6≤1.8 5.1.3 Degrees of freedom The four-bar mechanism is a bidimensional device so, in order to calculate the number of degrees of freedom (DOFs) the next equation must be used [32]: 𝑁º 𝐷𝑂𝐹′𝑠=3∗(𝑁−1)−2∗𝑃𝐼𝐼−𝑃𝐼 (5.2) Being N the number of bodies composing the mechanism (including the fix bar), 𝑃𝐼 the number of class-1- restrictions (those which remove one degree of freedom, for example, rotational pairs) and 𝑃𝐼𝐼 the number of class-2-restrictions (those which cancel two degrees of freedom, for example, prismatic pairs). It is direct to know that the 4-bar mechanism contains one DOF: 𝑁º 𝐷𝑂𝐹′𝑠=3∗(4−1)−2∗4=1 But in OpenSim, every mechanism or structure is defined in the tridimensional space, seven plane mechanisms as the 4-bar is considered to have 6 coordinates (3 translations and 3 rotations) instead of only 3 (2 translations and one rotation). This way, the equation to obtain the number of DOFs is the following one [32]: Inverse Dynamics in the 4-bar linkage 44 44 𝑁º 𝐷𝑂𝐹′𝑠=6∗(𝑁−1)−5∗𝑃𝑉−4∗𝑃𝐼𝑉−3∗𝑃𝐼𝐼𝐼−2∗𝑃𝐼𝐼−𝑃𝐼 (5.3) Being 𝑃𝐼, 𝑃𝐼𝐼, 𝑃𝐼𝐼𝐼, 𝑃𝐼𝑉 and 𝑃𝑉 those kinematic pairs which remove 1, 2, 3, 4 and 5 coordinates, respectively. If the number of DOFs is calculated for a 4-bar mechanism in the 3D space, this result is obtained: 𝑁º 𝐷𝑂𝐹′𝑠=6∗(4−1)−5∗4=−2 This negative result has no physical sense; it appears because there are redundant constraints in the mechanism, as OpenSim, as just explained, always considers the model to be tridimensional, assigning thus 6 coordinates to each body. This circumstance makes the model susceptible to lose its capacity to move due to numerical errors from software. That is the reason why it is needed to change any of the introduced pin joints for a ball joint, so as to cancel these redundant restrictions. In this project, the ball joint has been placed in the joint linking the input bar to the fix bar. 5.1.4 Generalized coordinates Given all introduced joints, finally 6 generalized coordinates define the movement of the 4-bar mechanism: • Θ12x, Θ12y and Θ12z, corresponding to the three relative spins between the input bar and the fix bar. In practice, as the motion is plane, Θ12x and Θ12y result to be negligible. • Θ23z, it is said, the relative rotation between the input bar and the coupler bar. • Θ34z, the relative spin between the coupler bar and the output bar. • Θ45z, the relative rotation between the coupler bar and the fix bar. It is also negligible, so the real number of generalized coordinates turns out to be 3. 5.2 Description of the motion As the 4-bar mechanism has only one degree of freedom, it is enough to use one generalized coordinate to express the movement of the model. In this case, as usual in problems involving this linkage, the coordinate related to the input angle, Θ12z, (hereafter, Θ12, as Θ12x and Θ12y are negligible) has been the chosen one. A 360 degrees spin has been imposed to this coordinate, parting from the following initial value: 𝜃2(0)=𝛱 𝑟𝑎𝑑 According to the definition of 𝜃2 given in Figure below, and corresponding to a vertical orientation of the input bar, which moreover spins, during a time span of 6 seconds, at constant velocity: 𝜔2=−1.0472 𝑟𝑎𝑑/𝑠 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 45 The inclusion of the negative sign means that the positive spin is counterclockwise, while the bar is spinning clockwise. This way, two phases can be distinguished in the motion: a) A first one, starting with the input bar in vertical position, with the initial condition right presented (it can be seen in Figure 5.2) where the bar is descending during its spin, aligning itself with the ground at the time instant t = 1.5 s (see Figure 5.3) until arriving at a new vertical position at t = 3 s (Figure 5.4). b) During the second phase of the motion, the bar performs the opposite motion, ascending in its rotation, aligning itself with ground at t =4.5 s (Figure 5.5) and turning to its initial position at t = 6 s (Figure 5.6). Figures 5.2 and 5.3. The 4-bar at t = 0 and at t = 1.5 s Figures 5.4 and 5.5. The 4-bar at t = 3 and at t = 4.5 s Figure 5.6. The 4-bar at t = 6 s Inverse Dynamics in the 4-bar linkage 46 46 5.3 Analytical inverse dynamics problem in the 4-bar mechanism The four bars mechanism is one of the paradigmatic examples in what refers to study multibody systems. It contains, as detailed before, four bodies linked by rotational joints. The Figure shows a scheme showing the defined generalized coordinates to express the motion: Figure 5.7. 4-bar mechanism built in WinMecC with generalized coordinates indicated [34]. Three generalized coordinates have been defined for each bar (without including the fixed bar): two describing the position of the center of gravity (ri = (xi yi )) and one for the orientation (θi). Thus, there are nine generalized coordinates in the model, ordered as shown in the expression below: 𝒒=(𝑥2 𝑦2 𝜃2 𝑥3 𝑦3 𝜃3 𝑥4 𝑦4 𝜃4) This manner, to calculate them, nine constraints must be defined, eight of them being geometrical constraints: 𝑥2−(𝐿2 2)∗cos(𝜃2)=0 (5.4) 𝑦2−(𝐿2 2)∗sin(𝜃2)=0 (5.5) 𝑥2+(𝐿2 2)∗cos(𝜃2)−𝑥3+(𝐿3 2)∗cos(𝜃3)=0 (5.6) 𝑦2+(𝐿2 2)∗sin(𝜃2)−𝑦3+(𝐿3 2)∗sin(𝜃3)=0 (5.7) 𝑥3+(𝐿3 2)∗cos(𝜃3)−𝑥4+(𝐿4 2)∗cos(𝜃4)=0 (5.8) Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 47 𝑦3+(𝐿3 2)∗sin(𝜃3)−𝑦4+(𝐿4 2)∗sin(𝜃4)=0 (5.9) 𝑥4+(𝐿4 2)∗cos(𝜃4)−𝐿1=0 (5.10) 𝑦4+(𝐿4 2)∗sin(𝜃4)=0 (5.11) While the ninth constraint is the kinematical one, imposing a time restriction to the input angle (in this case, subject to constant velocity ω =-2Π/6 rad/s): 𝜃2−𝜔∗𝑡=𝛱/2 (5.12) Being L1, L2, L3 and L4 the respective bar lengths, already defined in paragraph 6.1: 𝐿1=1.2 𝑚 𝐿2=0.4 𝑚 𝐿3=1 𝑚 𝐿4=0.8 𝑚 As detailed before, the motion consists of a 360-degrees spin of the input angle. Now, as explained in the mathematical definition, every generalized coordinate can be obtained. In this problem, those which arouse more interest are the rotational coordinates (consult Annex to see their graphic representations). Now, generalized velocities will be achieved, to that, the jacobian matrix must be calculated: 𝐶𝑞= ( 1 0 (𝐿2 2)∗sin(𝜃2)0 0 0 0 0 0 0 1 −(𝐿2 2)∗cos(𝜃2)0 0 0 0 0 0 1 0 −(𝐿2 2)∗sin(𝜃2)−1 0 −(𝐿3 2)∗sin(𝜃3)0 0 0 0 1 (𝐿2 2)∗cos(𝜃2)0 1 (𝐿3 2)∗cos(𝜃3)0 0 0 0 0 0 1 0 −(𝐿3 2)∗sin(𝜃3)−1 0 −(𝐿4 2)∗sin(𝜃4) 0 0 0 0 1 (𝐿3 2)∗cos(𝜃3)0 −1 (𝐿4 2)∗cos(𝜃4) 0 0 0 0 0 0 1 0 −(𝐿4 2)∗sin(𝜃4) 0 0 0 0 0 0 0 1 (𝐿4 2)∗cos(𝜃4) 0 0 1 0 0 0 0 0 0 ) (5.13) Also, the vector containing the derivatives with respect to time of the kinematical constraints is needed: 𝐶𝑡=(0 0 0 0 0 0 0 0 𝜔) (5.14) Inverse Dynamics in the 4-bar linkage 48 48 Now, using the Equation 4.3, velocities are obtained. As done with positions, rotational velocities result more important in this research (consult Annex to see their graphics). Now, to calculate accelerations, the first derivative respect to time of the jacobian matrix is needed: 𝐶𝑞󰇗= ( 0 0 (𝐿2 2)∗cos(𝜃2)∗𝜃2󰇗0 0 0 0 0 0 0 0 (𝐿2 2)∗sin(𝜃2)∗𝜃2󰇗0 0 0 0 0 0 0 0 −(𝐿2 2)∗cos(𝜃2)∗𝜃2󰇗0 0 −(𝐿3 2)∗cos(𝜃3)∗𝜃3󰇗0 0 0 0 0 −(𝐿2 2)∗sin(𝜃2)∗𝜃2󰇗0 0 −(𝐿3 2)∗sin(𝜃3)∗𝜃3󰇗0 0 0 0 0 0 0 0 −(𝐿3 2)∗cos(𝜃3)∗𝜃3󰇗0 0 −(𝐿4 2)∗cos(𝜃4)∗𝜃4󰇗 0 0 0 0 0 −(𝐿3 2)∗sin(𝜃3)∗𝜃3󰇗0 0 −(𝐿4 2)∗sin(𝜃4)∗𝜃4󰇗 0 0 0 0 0 0 0 0 −(𝐿4 2)∗cos(𝜃4)∗𝜃4󰇗 0 0 0 0 0 0 0 0 −(𝐿4 2)∗sin(𝜃4)∗𝜃4󰇗 0 0 0 0 0 0 0 0 0 ) (5.15) The rotation velocity is constant, so the term associated to the second derivative respect to time of the kinematical constraints is null. This way, the equation 4.5 can be applicated to obtain the accelerations. Those which correspond to rotational accelerations are displayed in Annex. The kinematic problem has been now solved for the 4-bar mechanism, so the dynamic problem can be solved. Starting with the mass matrix, the mass and the moment of inertia of each body must be defined. Every element has been assigned the same mass: 𝑚1=𝑚2=𝑚3=𝑚4=12 𝑘𝑔 With respect to the moment of inertia, the inertia matrix corresponding to an ellipsoid is the following one: 𝐼=1 5∗𝑚∗(𝑏2+𝑐20 0 0 𝑎2+𝑐20 0 0 𝑎2+𝑏2) (5.16) Being b the longest semiaxis, and a and c the other semiaxis. In this project, b is equal to each body’s length (already defined in subchapter 4.2), while a and c are equal to 0.01 m each. As the motion is represented in the OXY plane, only the third term of the inertia matrix results to be of interest. The Table 5.1 contain its value for each body in the four-bar mechanism: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 49 Table 5.1. Moment of inertia of each bar Body Input bar Coupler bar Output bar Izz (kg*m2) 0.08 0.05002 0.03202 Last, the vector containing forces must be defined. In this dynamic problem, only gravitational force was considered, so the vector would be the following one: 𝐹=𝑄𝑎𝑝𝑙+𝑄𝑣=(0 −𝑚2∗𝑔 0 0 −𝑚3∗𝑔 0 0 −𝑚4∗𝑔 0) (5.17) Where 𝑚2, 𝑚3 and 𝑚4 are the respective masses of bodies 2, 3 and 4, and g is the gravity, set to its value on Earth: 9,81 m/s2. Now, Lagrange multipliers may be obtained according to the equation 4.9. As the mechanism has one degree of freedom, in this case, the input angle θ2, the objective of the problem is to calculate the motor torque needed to perform the motion. This torque, named M2 as it is referred to the coordinate θ2, is associated to the kinematic constraint (the ninth one defined before) so, after applying the Equation 4.10 to obtain generalized forces and moments, the motor torque needed is the one corresponding to the last generalized force (associated to the ninth Lagrange’s multiplier): M2=−𝐶𝑞𝑇∗𝜆9 (5.18) Its graphic representation corresponds to the one displayed in Annex. The Inverse Dynamics problem is terminated now. The interest now is to calculate the dynamical power generated by the input bar, which is the only one that generates power, for being the only actuator in the model. The Equation serves to obtain this power: 𝑃=𝑀∗𝜔 (5.19) Being M the obtained motor torque and ω the angular velocity, described before. This result will be shown and commented in the next chapter. Once the analytical Inverse Dynamics problem applied to the 4-bar mechanism has been presented, it is the time to approach the problem from the point of view of OpenSim Moco. 5.4 MocoInverse: 4-bar Inverse Dynamics problem in Moco MocoInverse tool was presented in Chapter 3 as a way to perform Inverse Dynamics problems by solving an Inverse Dynamics in the 4-bar linkage 56 56 Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 57 6 RESULTS In this chapter, results which confirm hypotheses to prove in this work will be shown. To arrive to them, several intermediate results were necessary; these appear in the Annex. 6.1 Analytical Inverse Dynamics The first problem to solve was the analytical Inverse Dynamics in the non-muscled 4-bar, introduced in this work in paragraph 5.6. As explained then, the objective was to calculate the motor torque required to perform the motion, associated to the only degree of freedom in the mechanism, the input angle θ2. The result is displayed right below: Figure 6.1. Analytical motor torque of 4-bar. Obtained in Matlab. Regarding the power, this was the result calculated: Figure 6.2. Analytical power of 4-bar. Obtained in Matlab Results 58 58 Being this result of great importance in this work as it has served as the reference to compare the diverse values of power obtained with, so as to check the validity of methods to perform Inverse Dynamics. Observing the curve in Figure 6.2, and as expected knowing the two phases in which the motion is divided, also two phases can be clearly appreciated: a first one, spanning until the half of the cycle, where the power is negative, and a second one, parting from the time instant t = 3 s, where the power is positive. In the first half, the input is rotating downwards, so its objective is to oppose this fall, and exerts a negative power with this purpose. In the second phase, the mechanism is at its slowest position, and it needs to generate a positive power to raise the bar and turn it to the initial place. 6.2 Impossibility to use ID in closed-loop mechanisms in OpenSim GUI The first circumstance required to check was the impossibility to obtain correct results from the Inverse Dynamics performed in the Graphical User Interface of Moco in a closed-loop mechanism. This can be checked by comparing the dynamical equivalent power obtained in this problem to the power resulting in the analytical problem. These graphics are shown in Figures 6.3 and 6.4, respectively: Figure 6.3. Power resulting in ID by OpenSim. Obtained in Matlab. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 59 Figure 6.4. Analytical power of 4-bar. Obtained in Matlab It is clearly seen that these powers do not coincide, what verificates that OpenSim GUI is not a valid tool to perform Inverse Dynamics. This circumstance was waited, because the reason why OpenSim cannot solve Inverse Dynamics in closed-loop mechanisms is that it does not recognize constraints in the dynamical problem, what leads to errors at the moment of calculating forces and moments. 6.3 Capacity of OpenSim Moco to solve ID in closed-loop mechanisms The next step consisted of determinate if the problem solved by using MocoInverse, in the closed-loop mechanism (in this case, the 4-bar without muscles) offered correct and reliable results, it is said, the power generated in this problem might be the same as in the analytical one. Figure 6.5 displays the power generated by the 4-bar mechanism without muscles, and Figure 6.6 compares it to the obtained in the analytical case: Figure 6.5. Power in a non-muscled 4-bar mechanism. Obtained in Matlab. Results 60 60 Figure 6.6. Comparison of total power obtained in Moco and the analytical way. Obtained in Matlab. Results confirm what was being tried to be proven: the power generated by the 4-bar is the same calculated by the analytical way and by using Moco. This confirms the validity of this tool to perform Inverse Dynamics analysis in closed-chain mechanisms. 6.4 Capacity of OpenSim Moco to solve ID in muscled closed-loop mechanisms Now it is the moment to evaluate if OpenSim Moco can serve to study closed-loop mechanisms with muscles, assimilable to models involving human people wearing exoskeletons. The requirement is still that the power supplied by the muscled 4-bar were the same as in the model without muscles, because, as explained before, muscles play the role of actuators, providing force (and thus, torque and power), but the motion is still the same: the input bar must complete a 360-degrees rotation, and the analytical solution informs about the required torque to fulfill this movement. Figure 6.7. Power resulted in muscled 4-bar. Obtained in Matlab. Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 61 Figure 6.8. Comparison of total power obtained in Moco and the analytical way, in the muscled 4-bar. Obtained in Matlab. In Figure 6.7, power supplied by the muscled 4-bar appears, while in Figure 6.8 it is compared to the analytical one. It is clear that, effectively, the power calculated in the problem involving the muscled 4-bar mechanism is the same as the one obtained in the analytical problem. Therefore, OpenSim has been confirmed as an accurate tool to exert Inverse Dynamics problems in closed-loop musculoskeletal models, in which problems involving exoskeletons are, as previously said in this work, included. But analysis of results must continue yet. Now, the influence of muscles in the motion will be detailed, focusing on their activations, the forces they generate, both passive and active, and how they contribute to the total power in contrast to the Coordinate Actuators. In first place, the contribution of these Coordinate Actuators and the muscles to the total power is analyzed in Figure 6.9: Figure 6.9. Contribution of muscles and Coord. Actuators to power. Obtained in Matlab. Results 62 62 It is observed that the muscles (whose power is blacklined) cannot generate the whole needed power to perform the motion. In some parts of the cycle, they exert the required part (from t = 0 to t = 0.5 s and from t = 3 s to t = 5.7 s. During the rest of time, Coordinate Actuators have to supply power too, and even in two instants (t = 1.5 s and t = 4.5 s) these actuators provide the 100% of the power. The explanation to that resides on the fact that the input bar is aligned to the ground in those instants, so the moment arm to generate the torque is zero. Muscles cannot therefore produce torque nor power. The appearance of the Coordinate Actuators in the motion is directly related to muscle activation. It can be observed in Figure, where powers have been normalized so as to focus on the evolution of the graphic and not in magnitude orders: Figure 6.10. Comparison of powers by muscle and their activations. Obtained in Matlab It is observed that muscles are able to produce the whole force required until they reach their maximum level of activation. This means they cannot generate more force, and so, the Coordinate Actuators start to supply torque. When muscle activation is under 1, muscles are generating all the force needed, and they do not need to get more activated. In order to study the relation of the activation with respect to the active and passive forces, as well as the muscle parameters intervening in the force (fiber length and contraction velocity) the next Figures will be useful: Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 63 Figure 6.11. Right muscle forces. Obtained in Matlab Figure 6.12. Left muscle forces. Obtained in Matlab It can be observed in both muscles that passive forces appear when the fiber lengths surpass its optimal fiber length. This can happen even if the muscle is not activated, for example, as it happens in the left muscle between the instants t =1.5 s and t = 3.3 s. The muscle is shortening, but it will be producing passive forces while its fiber length is over its optimal length. On the other hand, there is no active force if the activation is zero, as expected. But the value of the active force also depends on the active force-length-velocity parameters, as it happens in right muscle between instants t =4 s and t=4.5 s, when the active force grows up since these parameters are growing up. When the fiber length gets closer to its optimal length, higher is the value of the active force-length parameter, the same occurs with force velocity parameter when the contraction velocity gets increased. But when fiber length is too much bigger or smaller than its optimal length, or contraction velocity is too short (or even negative, what means that muscle is lengthening), the active force length gets decreased and so, the muscle can generate less power. Results 64 64 To sum up, it can be said that muscles stop generating the force required by the problem when their activation reaches its maximum value, and consequently, Coordinate Actuators supply power. Furthermore, if the fiber length is far from its optimal value, and contraction velocity is small (or negative), the production of force by the muscle will be limited. This situation can be improved if the position of the muscle is changed, what will be seen in the next paragraph. 6.5 Influence on power of changes in the insertion point of muscles Last, to close this work, a sensitivity analysis has been done to study what happens to muscle power if the insertion point is changed. This way, the insertion points of muscles in the input bar have been modified the way shown in Figures 6.13 and 6.14: Figures 6.13 and 6.14. New insertion points of muscles in input bar. OpenSim [6] In the first case, the insertion point common to both muscles has been placed above the center of gravity of the input bar; concretely, in the middle point between the center of gravity and the vertix common to the input and the coupler bar; in the second case, the insertion point has been collocated between the center of gravity and the vertix of the input bar in contact to the ground. Then, in the third phase of this analysis of sensitivity, the two muscles inserted in the center of gravity of the input bar have been maintained, and the two muscles situated above have been have been considered, as shown in Figure 6.15: Figure 6.15. 4-muscles 4-bar linkage in OpenSim [6] Results of this analysis of sensitivity are now shown. Starting with the study of influence of the collocation of Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 65 the insertion point, the graphics corresponding to the power generated by muscles is displayed in Figure 6.16, while the one performed by the actuators can be observed in Figure 6.17: Figure 6.16. Influence of insertion point in muscle power. Obtained by Matlab. Figure 6.17. Influence of insertion point in actuators’ power. Obtained by Matlab. It can be seen that locating the insertion point above the middle point of the input bar (this case has been represented by the green line) results favorable during most of the cycle in terms of resemblance of muscle power to total power, because the power generated by the muscles equals the power required to exert the motion during more time than in the case where the insertion point is in the center of gravity of the bar. As the power generated by Coordinate Actuators is the resulting from subtract muscle power from total, the power from actuators is zero during these instants of time. This result can seem logical, because at placing the insertion point this way, the moment arm gets increased, so it is easier for the muscle to generate more torque. Bibliography 72 72 [26] Nagano, A., Gerritsen, K.G.M.: Effects of neuromuscular strength training on vertical jumping performance—a computer simulation study. J. Appl. Biomech. 17, 113–128 (2001). [27] OpenSim, OpenSim Documentation. First Order Activation Dynamics. October 2021. Available online: https://simtk-confluence.stanford.edu:8443/display/OpenSim/First-Order+Activation+Dynamics [28] OpenSim, OpenSim Documentation. Characteristic Musculotendon curves. October 2021. Available online: https://simtk-confluence.stanford.edu:8443/display/OpenSim/Characteristic+Musculotendon+Curves [29] OpenSim, OpenSim Documentation. How Inverse Kinematics works. October 2021. Available online: https://simtk-confluence.stanford.edu:8443/display/OpenSim/How+Inverse+Kinematics+Works [30] OpenSim, OpenSim Documentation. How Inverse Dynamics works. October 2021. Available online: https://simtk-confluence.stanford.edu:8443/display/OpenSim/How+Inverse+Dynamics+Works [31]. Exoskeleton model in OpenSim. Departamento de Ingeniería Mecánica y Fabricación. Universidad de Sevilla. 2020. [32]. Chamorro Moreno, R. et al., Teoría de Máquinas y Mecanismos. Departamento de Ingeniería Mecánica y Fabricación. Universidad de Sevilla. 2015. [33]. Escalona Franco, J.L., Cinemática y Dinámica de Máquinas. . Departamento de Ingeniería Mecánica y Fabricación. Universidad de Sevilla. 2015. [34]. WinMecc. Área de Ingeniería Mecánica de la Universidad de Málaga. [35] OpenSim Moco, OpenSim::MocoControlGoal. October 2021. Available online: https://opensimorg.github.io/opensim-moco-site/docs/ [36] OpenSim Moco, OpenSim::MocoSunSquaredStateGoal. October 2021. Available online: https://opensimorg.github.io/opensim-moco-site/docs/ [37] OpenSim Moco, OpenSim::MocoInitialActivationGoal. October 2021. Available online: https://opensimorg.github.io/opensim-moco-site/docs/ Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 73 Annex. Intermediate results 74 74 ANNEX. INTERMEDIATE RESULTS In this annex, some intermediate results without direct relationship with the purpose of this work are mostrated. They are several muscle parameters obtained during the process to calculate passive and active forces (described in Chapter 4). Figure 9.1. Normalized tendon force in both muscles. OpenSim Moco [13] Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 75 Figure 9.2. Normalized tendon force in both muscles. OpenSim Moco [13] Figure 9.3. Normalized tendon length in right muscle. OpenSim Moco [13] Annex. Intermediate results 76 76 Figure 9.4. Normalized tendon force in left muscle. OpenSim Moco [13] Figure 9.5. Muscle lengths. OpenSim Moco [13] Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 77 Figure 9.6. Muscle fiber lengths. OpenSim Moco [13] Figure 9.7. Normalized contraction velocities in both muscles. OpenSim Moco [13] Annex. Intermediate results 78 78 Figure 9.8. Normalized passive force in both models, right muscle. OpenSim Moco [13] Figure 9.9. Normalized left passive force in both models. OpenSim Moco [13] Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 79 Figure 9.10. Right active force-length multiplier in both models. OpenSim Moco [13] Figure 9.11. Left active force-length multiplier in both models. OpenSim Moco [13] Annex. Intermediate results 80 80 Figure 9.12. Right active force-velocity multiplier in both models. OpenSim Moco [13] Figure 9.13. Left active force-velocity multiplier in both models. OpenSim Moco [13] Analysis of feasibility in OpenSim Moco of Inverse Dynamics problems involving human models with exoskeletons 81 Figure 9.14. Normalized active forces in both muscles. OpenSim Moco [13]