Development of muscle torque generators for optimal control tracking of human walking

Treball de Fi de Grau Grau en Enginyeria en Tecnologies Industrials (GETI) Development of muscle torque generators for optimal control tracking of human walking MEMÒRIA Autor: David Lasierra Hernández Directors: Míriam Febrer Nafría Josep Maria Font Llagunes Convocatòria: Juliol 2022 Escola Tècnica Superior d’Enginyeria Industrial de Barcelona Development of MTGs for OCT of human walking Page 1 Abstract Muscle torque generators (MTGs) can be used as an alternative of detailed musculoskeletal models in computer simulations of human movement, and have been used by researchers to generate accurate dynamic simulations while mantaining a reasonable degree of biofidelity. The main goal of this thesis is to improve the reliability of two torque-driven models (i.e., at skeletal level) of different complexity when muscle torque generators are added to them in order to perform a concrete motion. For this purpose, we defined two problems based on different optimal control techniques. The first problem predicts the motion of a simple pendulum model and the second one tracks experimental data from a human walking motion in a 2D HAT (head, arms and trunk) model. In both problems, two versions of the models are considered: a torquedriven version and a version in which muscle torque generators are employed. Both torquedriven models are obtained from the OpenSim repository and the implementation of the MTGs has been programmed with MATLAB. For each problem, different studies were conducted to reach the stated goal. In the predictive problem, several simulations were carried out based on different initial and final conditions. Subsequently, a qualitative evaluation between both models showed that muscle torque generators were slightly advantageous than its torque-driven version when the studied simulation encompassed a shorter range of motion. Conversely, in the tracking problem, a more complex study was performed, based on different optimal control formulations. A final quantitative assessment showed that the tracking in torques was adequate and quite similar between models, but when tracking experimental coordinates the torque-driven version provided better results. Despite the fact the tracking with MTGs was correct, a more realistic behaviour than its torquedriven version was expected. This project can be considered a first research work in muscle torque generators modelling in the Biomechanical Engineering Lab (BIOMEC) at UPC. Thus, further studies would need to be conducted in order to obtain a more realiable modelling. Parameters such as the maximum isometric torque or the characteristic musculotendon curves involved in the MTGs could be more in depth investigated to gain better insights. Development of MTGs for OCT of human walking Page 2 Contents 1 Introduction 6 1.1 Motivation ......................................... 6 1.2 Project objectives ..................................... 6 2 Theoretical background 7 2.1 Biomechanics of human motion ............................. 7 2.1.1 Anatomical planes ................................ 7 2.1.2 Gait cycle ..................................... 8 2.2 Multibody system modelling .............................. 9 2.2.1 Skeletal modelling ................................ 9 2.2.2 Individual muscle modelling .......................... 10 2.2.3 Muscle torque generators ............................ 11 2.2.4 Foot-ground contact modelling ......................... 13 2.3 Dynamics Analysis .................................... 13 2.3.1 Inverse Dynamic Analysis ............................ 13 2.3.2 Forward Dynamic Analysis ........................... 14 2.4 Optimal control prediction ................................ 15 3 Methods 16 3.1 Simple pendulum simulations .............................. 16 3.1.1 MTGs modelling ................................. 17 3.1.2 Optimal control predictive problem ...................... 18 3.1.3 Simulations of interest .............................. 19 3.2 2D HAT simulations ................................... 20 3.2.1 MTGs modelling ................................. 21 3.2.2 Optimal control tracking problem ....................... 24 3.2.3 Simulations of interest .............................. 24 4 Results and discussion 26 4.1 Simple pendulum predictive problem ......................... 26 4.1.1 Simulation I .................................... 26 4.1.2 Simulation II ................................... 27 4.1.3 Simulation III ................................... 28 4.1.4 Simulation IV ................................... 28 4.1.5 General discussion ................................ 29 4.2 2D HAT tracking problem ................................ 30 4.2.1 Study of different formulations in the MTG cost function .......... 30 4.2.2 Influence of passive elements and joint damping in motion tracking . . . . 31 4.2.3 Evaluation of muscle torque generators implementation .......... 31 4.2.4 General discussion ................................ 34 5 Project impact 35 5.1 Economic cost ....................................... 35 5.2 Environmental impact .................................. 36 Conclusions 37 Acknowledgements 38 3 List of Figures 1 The anatomical planes from the human body. Extracted from [17]. ........ 7 2 Human walking gait cycle of the right leg. Extracted from [26]. .......... 8 3 Stride describing the gait cycle. Extracted from [26]. ................. 8 4 2D biomechanical model used in a running gait analysis. Extracted from [25]. . . 9 5 The Hill-type muscle model. Adapted from [9] and [37]. .............. 10 6 The characteristic musculotendon curves. Extracted from [9]. ........... 11 7 Biceps and triceps as an example of antagonistic pairs. Adapted from [3]. . . . . 11 8 Model actuated by (a) Hill-type muscles and (b) MTGs. Extracted from [15]. . . 12 9 Foot-ground contact modelled as (a) a kinematic constraint between an ellipse and a plane (extracted from [18]) and (b) a constitutive volumetric contact model. 13 10 Diagram of an inverse dynamic analysis. ........................ 14 11 Diagram of a forward dynamic analysis. ........................ 14 12 Simple pendulum configuration as a (a) torque-driven model (b) MTGs-based model. ........................................... 16 13 Simple pendulum simulations. ............................. 19 14 Bodies and generalized coordinates of the model. .................. 20 15 2D HAT configuration as a (a) torque-driven model (b) MTGs-based model. . . 21 16 Process that leads to the final version of the MTG model. .............. 25 17 Evolution of states and joint torque in the torque-driven model and the MTGdriven model when Simulation I is conducted. .................... 26 18 Evolution of states and joint torque in the torque-driven model and the MTGdriven model when Simulation II is conducted. .................... 27 19 Evolution of states and joint torque in the torque-driven model and the MTGdriven model when Simulation III is conducted. ................... 28 20 Evolution of states and joint torque in the torque-driven model and the MTGdriven model when Simulation IV is conducted. ................... 29 21 Tracking of coordinates and torques in a 2D HAT model based in a torque-driven (in blue) and an MTG (in red) versions. Experimental data is shown with a dashed black line. The gait accounts the 80% of the cycle and left joints are shown. 32 22 Pairs of muscle torque generators at each left joint and lumbar. ........... 33 List of Tables 1 Simple pendulum optimal control predictive statement. ............... 18 2 Coordinates and DoF. ................................... 20 3 Parameters that define the characteristic curve of the MTG in the 2D HAT problem. 23 4 2D HAT Torque-driven optimal control statement. .................. 24 5 2D HAT MTG-driven optimal control statement. ................... 24 6 Set of cost functions studied. ............................... 25 7 RMSE of coordinates obtained from the different formulations. ........... 30 8 RMSE of torques obtained from the different formulations. ............. 30 9 RMSE of coordinates with and without passive elements (PE) and joint damping (DAMP). .......................................... 31 10 RMSE of torques with and without passive elements (PE) and joint damping (DAMP). .......................................... 31 11 RMSE of coordinates in the torque-driven model (TD) and the MTG model. . . . 31 12 RMSE of torques in the torque-driven model (TD) and the MTG model. ..... 32 4 Development of MTGs for OCT of human walking Page 5 13 Computational time and number of iterations required in both models ...... 34 14 Calculation of the final project cost. Variable costs of the computer and licenses are obtained from dividing their fixed cost by their useful life in hours. Variable cost of electrical energy is found from its price multiplied by the power consumption of computers and light. ............................... 35 Development of MTGs for OCT of human walking Page 6 1 Introduction 1.1 Motivation The human musculoskeletal system, responsible for the locomotion, is composed by about 200 bones and 300 muscles. The development of models of such a complex system is a non-trivial task that involves the traditional compromise between straightforwardness and accuracy. On the one hand, the models should be complex enough to deliver accurate estimations of the phenomena studied. On the other hand, the models should be simple enough to keep the problem tractable and the results and analysis straightforward [1]. Using detailed musculoskeletal models in computer simulations of human movement can provide insights into individual muscle and joint loading; however, these muscle models increase problem dimensionality and require difficult-to-fit parameters [15]. This prompts the question regarding what methods are available to maintain accuracy within a given musculoskeletal model while limiting the required computational demand. To mitigate these effects while maintaining a reasonable degree of biofidelity, we employ muscle torque generators (MTGs). The implementation of MTGs in complex models can be tedious at first. However, if the model is simplified from several degrees of freedom to only one, then the problem formulation becomes easier to understand and also an extension of the problem to a more complex version can be easily done. Therefore, in this thesis, two problems of different complexity have been evaluated. In both cases, a torque-driven model and a model which includes MTGs have been studied. 1.2 Project objectives The main goal of this thesis is to improve the reliability of two torque-driven models of different complexity when muscle torque generators are added to them in order to perform a concrete motion. The first study is based on a predictive optimal control problem that uses a simple pendulum as a model. Whereas the second study is based on an optimal control tracking problem, in which the ability of the model to reliably reproduce a gait cycle is evaluated. The specific objectives that involve the main goal of this project are: •Research the state of the art in the modelling of muscle torque generators for human motion simulations. •Implement muscle torque generators in a simple pendulum model with MATLAB. •Compare the predictive behaviour between the torque-driven model and the MTG-based model through a qualitative evaluation. •Study different optimal control formulations for human gait tracking when muscle torque generators are taken into account in the model. •Implement muscle torque generators in a 2D HAT model with MATLAB. •Compare the tracking behaviour between the torque-driven model and the MTG-based model when a quantitative evaluation is conducted. Development of MTGs for OCT of human walking Page 7 2 Theoretical background This section gives theoretical insight into the basis of human walking and its modelling, as well as some key concepts used throughout the project. 2.1 Biomechanics of human motion 2.1.1 Anatomical planes In order to describe the orientation and location of human body structures it is convenient to define some reference planes (see Figure 1). They are called anatomical planes and they separate the body into different sections. Those are: •The sagittal plane. It splits the human body up into left and right sections. •The coronal plane. It divides the body into posterior and anterior portions. •The transverse plane. It separates the body into an upper and a lower part. Figure 1: The anatomical planes from the human body. Extracted from [17]. Depending on the type of motion that is analyzed, some planes contain more relevant information than others. Even so, many textbook authors [35] and researchers emphasize in the use of the sagittal plane in gait analysis and ignore the other two. However, when a motion from a pathological subject is analyzed, other planes (e.g. the coronal plane in the case of bilateral hip pain) would yield relevant information as well [5]. In this project, a 2D gait model will be used. Only motion in the sagittal plane, which is where most of the movement takes place, will be considered. Development of MTGs for OCT of human walking Page 14 Mathematically speaking, this approach requires the system’s coordinates and its derivatives. Therefore, by evaluating them in the motion equations for a time range, joint forces and torques can be computed. Reflective markers and force plates are used to obtain motion data and GRF respectively. Since these devices do not provide directly the coordinates nor their time derivatives, a preprocess called inverse kinematics (IK) must be applied (see Figure 10). This consists in solving an optimization problem where the difference between captured marker trajectories and model marker trajectories is minimized. Figure 10: Diagram of an inverse dynamic analysis. 2.3.2 Forward Dynamic Analysis Forward dynamic analysis describes the motion of the system when some particular forces or torques are applied. The input data in this approach are the muscle forces or the resultant joint torques. Body segment parameters and GRF are considered as inputs as well (see Figure 11). It is also possible to find muscle excitations as input data. In that case, the simulation is called muscle-driven simulation. Unlike the IDA, in FDA the differential equations of motion shall be integrated with respect to time in order to find the evolution of joint coordinates. This can sometimes lead to integration errors and thus, combined with the unstable character of human walking, conduct to non-stable solutions. However, these errors can be compensated with the use of control methodologies. Figure 11: Diagram of a forward dynamic analysis. In this thesis, an inverse dynamic analysis is carried out within an optimal control formulation. This formulation is explained in the following section 2.4. Development of MTGs for OCT of human walking Page 15 2.4 Optimal control prediction Techniques defined in the previous section 2.3 can be used to predict novel human motions (a priori unknown) for which no experimental data are available at first. During the last years, there has been a growing interest in motion prediction due to its applications. Some of them are: anticipating the result of a surgery, designing assistive devices or analysing the dynamic simulation of a specific motion [24]. There are several ways to predict human motion. However, the most accepted approach consists of formulating an optimisation problem. The main idea is to find which system configuration provides an optimal solution to a problem given some constraints and a cost function to be optimized. Therefore, any optimal control problem is usually defined by the following parameters: •State variables. They are used to describe the current state (position, velocity) of the system. Examples of state variables in gait simulations are: joint coordinates, velocities, muscle length... •Control variables. These variables manipulate the state variables in order to satisfy some desired conditions. Examples of control variables in gait simulations are: joint accelerations, muscle activation/excitation, muscle torques... •Dynamic constraints. They are a set of equations governed in the states and have to be satisfied during the optimization problem. •Path constraints. They are inequalities that must be satisfied during an entire time interval. These restrictions can either be bounds on states and controls, or algebraic path constraints. •Boundary conditions on state and time. Some state and time situations are imposed at the beginning and at the end of the time interval. •Cost function. It is the mathematical expression that wants to be optimized, which is a function of state and control variables. Different cost functions have been used to predict gait. Some examples are minimizing weighted normalized torques [29] or joint accelerations as well as minimizing muscle activations [13]. In this thesis, we employ optimal control within an inverse dynamic analysis in the problem formulation. The two problems studied in this project (see sections 3.1.2 and 3.2.2) have a similar structure in terms of the problem formulation. However, the first one predicts a novel motion and the other tracks a motion from experimental data. The main difference in these two kinds of formulation lies in the cost function. Development of MTGs for OCT of human walking Page 16 3 Methods This section introduces the models used throughout the project as well as how the muscle torque generators have been implemented in each of them. The implementation of MTGs in complex models can be tedious at first. If the model is simplified from several degrees of freedom to only one, then the problem formulation becomes easier to understand. Not only is it easier to implement the MTGs, but also an extension of the problem to a more complex version can be easily done. In order to ease the understanding of MTGs in human walking simulations, a simpler problem has been defined. Thus, it was considered that the formulation of a simple pendulum with a single joint could be easily transferred to several joints of the human body. Hence, the chapter is divided into two different parts. The first part corresponds to an optimal control predictive problem using a simple pendulum model, while the second one is an optimal control tracking problem using a 2D HAT model. That means that in the simple pendulum problem we look for a new motion and in the 2D HAT problem we look for a motion which is as close as possible to an experimentally measured one. In each of the parts two different models are presented: a torque-driven model (i.e., at skeletal level) and a model where muscle torque generators have been implemented. In both problems, the torque-driven model is obtained from the OpenSim repository and the implementation of the MTGs has been programmed with MATLAB. 3.1 Simple pendulum simulations The simple pendulum problem consists of an m= 1 kg massive bob which is attached to a massless stick with length ℓ= 0.5mand they both swing back and forth in a periodic motion. The joint where the motion is produced is fixed to the ground and hence, only rotational motion is allowed. The coordinate that describes this one-degree-of-freedom problem is the angle qand it will change from an initial position qini to a final one qfin within a bounded time t∈[tini, tfin]. The motion of the pendulum in the torque-driven problem is given by an external moment Mz (see Figure 12a), while in the problem with MTGs, it is due to two torques (see Figure 12b), one associated with the extension τME and the other referred to flexion τMF . Figure 12: Simple pendulum configuration as a (a) torque-driven model (b) MTGs-based model. Development of MTGs for OCT of human walking Page 17 3.1.1 MTGs modelling The net torque at the model’s joint is the sum of the signed flexor and extensor muscle torques acting at that joint: τ=τMF +τME (1) Since each muscle torque acts in a single direction, there are two MTGs acting at the joint, a flexor τMF and an extensor τME. The torque τMdeveloped by a single MTG is given by an adapted formulation from [19]: τM=τM oafA(q)fV( ˙q)(2) The torque developed by Equation 2is a function of the control input afrom the solver, the angle q, and the angular velocity ˙qof the joint. While the angle of the joint changes the value of the passive torque-angle curve fA(q), the angular velocity of the joint affects the value of the torque velocity curve fV( ˙q). Finally, the parameter τohelps to build the characteristic curve of the MTG as it provides the maximum value achieved during the motion. Either the flexor and extensor maximum torques were set to τMF o=τME o= 10 Nm. These values were adjusted from the torque-driven problem. In this model, the characteristic curve fA(q)was defined as a parabolic function as it follows: fA(q) = 1 4p·q2−h 2p·q+h2+ 4pk 4p(3) where (h, k)correspond to the vertex coordinates and pchanges the concavity of the function. These parameters were set as: h=−0.02 rad,k= 1 and p=−0.57 rad2. Conversely, the torque velocity curve fV( ˙q)was extracted from the literature [28] and can be written as the following piecewise function: fV( ˙q) =                    0for ˙q≤ −1, 1+ ˙q 1−q kCE1 for −1≤˙q≤0, 1+ ˙qfmax kCE2 1+ ˙q kCE2 for ˙q > 0 (4) where kCE1and kCE2are force velocity shape factors and are set to kCE1= 0.25 and kCE2= 0.06 in this work. In contrast, fmax is the maximum normalized achievable torque and it is set to fmax = 1.6. Development of MTGs for OCT of human walking Page 18 3.1.2 Optimal control predictive problem To predict the simple pendulum motion, two different models have been used: the torquedriven model and the MTG-based model. Therefore, the formulation of the problem differs depending on which model is studied. Table 1shows the optimal control configuration for each of the models, where all the magnitudes are represented in the international system units. Both problems have the same state variables and the angular coordinate qwill move from an initial position qini to a final one qfin within a bounded time t∈[tini, tfin]. Note that a reserve actuator τreserve was added in the controls from the MTG-driven statement. This torque will help to satisfy the path constraints when the MTG torques are unable to perform the desired motion. Also, note that the activations are only considered in the MTG-driven statement since they are only dependant on the muscles. As it is stated in section 2.4, these activations together with Mzas well as τreserve will be considered in the cost function to minimise. It is also remarkable that the path constraints are different for each model. While in the torquedriven model the generated torque at the joint has to be equal to the one computed with the inverse dynamic analysis, in the MTG-driven is the sum of the flexor and extensor muscle torques (plus the reserve actuator). Table 1: Simple pendulum optimal control predictive statement. Torque-driven MTG-driven States q,˙q q,˙q State constraints −π 2≤q≤π 2−π 2≤q≤π 2 −100 ≤˙q≤100 −100 ≤˙q≤100 Controls ¨q, Mz¨q, τreserve, a Control constraints −100 ≤¨q≤100 −100 ≤¨q≤100 −10 ≤Mz≤10 −10 ≤τreserve ≤10 0≤a≤1 Path constraints τ=Mzτ=τMF +τME +τreserve Cost function J= 0.1τ2+ 0.01¨q2J= 0.1τ2+a2+τ2 reserve In order to implement the optimal control algorithm, we used GP OP S−II, which works within aMAT LAB environment. Development of MTGs for OCT of human walking Page 19 3.1.3 Simulations of interest For each of the two models, four different simulations have been carried out based on different initial and final conditions (see Figure 13). The main purpose of these has been twofold. On the one hand, we wanted to check the consistency of the MTGs for each simulation. For this purpose, the evolution of the system states (q, ˙q) as well as the Mztorque will be qualitatively compared in both models. On the other hand, a study will also be made comparing the four simulations in which the MTGs have been implemented. In this way, it will be possible to check if the pendulum shows any kind of tendency in its motion depending on its initial and final conditions when the MTGs are applied. Figure 13: Simple pendulum simulations. Note that each of the simulations has been conducted in both the torque-driven version and the version with muscle torque generators implemented. Development of MTGs for OCT of human walking Page 20 3.2 2D HAT simulations Having studied the simple case for a single joint as was the case of the pendulum, we can now proceed to a more complex study, which will include several joints as is the case of human walking. Bodies and generalized coordinates of the model are illustrated in Figure 14. The studied problem consists of ten degrees of freedom (see Table 2). From these ten, the three involving the pelvis are actuated by residual forces and moments. These have been artificially added so that the equations of motion of the system are dynamically consistent. Table 2: Coordinates and DoF. Coordinate Degree of freedom (DoF) x0Pelvis xdisplacement y0Pelvis ydisplacement q0Pelvis tilt q1Lumbar extension q2Right hip flexion angle q3Right knee angle q4Right ankle angle q5Left hip flexion angle q6Left knee angle q7Left ankle angle Figure 14: Bodies and generalized coordinates of the model. The approach that has been studied in this part of the project is to add muscle torque generators on the OpenSim existing model, which has been scaled to the subject for which the experimental data were available.The main objective is to study the consistency of the model when MTGs are added to it. For this purpose, two studies will be carried out: a torque-driven version (see Figure 15a) and a version where MTGs have been implemented (see Figure 15b). Similarly to the pendulum, the torque generators in the first model are external moments τiat each joint i∈ {1,2, ..., 7}; whereas in the second model a pair of muscle torque generators τMF iand τME iare acting at each joint. Note that both models contain the residual forces and moment acting in the pelvis, but neither an external moment nor a pair of MTGs are applied in it such as in the rest of joints. Development of MTGs for OCT of human walking Page 21 Figure 15: 2D HAT configuration as a (a) torque-driven model (b) MTGs-based model. Finally, both models will be compared with experimentally obtained data, which include the coordinates qexp as well as the torques τexp involved in the tracking. These data were obtained from markers and force plates respectively and were collected in a previous study [22]. Therefore, in this problem, an optimal control framework that tracks experimental data will be applied in both 2D HAT models. 3.2.1 MTGs modelling In this study, the human body is modelled as a sagittal-plane multibody system that is actuated by agonist and antagonist pairs of muscle torque generators at each joint. In comparison to the simple pendulum problem, apart from studying several joints, other parameters in relation to the suppression of vibrations in the model as well as the incorporation of passive elements will be considered. Then, the net torque at each of the model’s joints is the sum of the signed flexor and extensor muscle torques acting at that joint and joint damping [19]: τi=τMF i+τME i−β˙qi(5) Since each muscle torque acts in a single direction, there are two MTGs acting at the joint, a flexor τMF iand an extensor τME i, for a total of 14 MTGs for the whole model. The parameter βcorresponds to joint damping coefficient defined at Equation 7and ˙qiis the derivative of the coordinate qi(i.e., the angular velocity acting at the joint i). Development of MTGs for OCT of human walking Page 22 The torque τMdeveloped by a single MTG is given by [19]: τM=τM o afA(q)fV( ˙q) + fPE (q)1−βP E ˙q ˙qM max !(6) The torque developed by equation 6is a function of the control input afrom the solver (in this case mapped to the activation of the muscle), the angle q, and angular velocity ˙qof the joint. The angle of the joint changes the value of fA(q), the passive torque-angle curve. The angular velocity of the joint affects the value of fV( ˙q), the torque velocity curve, and also the damping torque of the passive element. A non-linear normalized damping term βPE is added to the passive element to suppress possible vibrations. This parameter is usually set to 0.1 in muscles and thus, it was also decided to keep that value in this model. Another parameter that is taken into account when passive elements are considered is ˙qM max, which corresponds to the maximum angular velocity achieved by the joint. Note that this value changes depending whether the torque is the flexor or the extensor. Finally, τM ois the maximum isometric torque achieved by the joint. Similarly to the maximum angluar velocity, τM ois also dependant on the nature of the torque and thus, it will be different for both the flexor and extensor torques. The light damping at the joint in Equation 5is the passive damping introduced by the musculature and tissue surrounding the joint. The damping coefficient is defined as [19]: β=ητMF o+τME o ˙qMF max + ˙qME max (7) Therefore, the amount of damping is proportional to the strength of the musculature and inversely proportional to its maximum angular velocity. The parameter ηis a normalized joint damping scaling factor and it can take the values 0.2 or 0.4 depending on which kind of joint is studied. The first value corresponds to lower body joints and the latter to arms. Since the model used in this thesis only accounts on lower body joints, then ηis set to 0.2. Several literature sources are used to build the characteristic curves for the MTGs (see "MTG Parameters" at Table 3). Regarding the musculotendon characteristic curves, the torque-angle and torque-velocity curves have been described as in the pendulum problem (see Equation 3 and Equation 4) but transferred to several joints. While the parameters p, h and kthat define the first mentioned curve change according to Table 3(refer now to "fParameters"), the parameters that define the torque velocity curve kCE1, kCE2and fmax are kept as in the pendulum formulation 3.1.1. Finally, the passive torque angle curve is defined according to the literature [28] as: fPE (q) = e(kP E ·q−l0)/εM 0−1 ekP E −1(8) The parameter kPE is a shape factor, εM 0is the parallel element strain and l0is the slack angle which is different for each joint. These parameters are represented at Table 3as well. Development of MTGs for OCT of human walking Page 23 Table 3: Parameters that define the characteristic curve of the MTG in the 2D HAT problem. MTG Parameters Lumbar Hip Knee Ankle References Extension Flexion Extension Flexion Extension Flexion Extension Flexion τ0(Nm) 275.1 211.7 175.7 157.3 285.6 98.60 127.6 44.30 Millard 2017 [19] ˙qmax (rad/s) 0.38 0.65 12.00 12.40 16.60 22.40 8.00 10.70 JacksonMI [16] βPE 0.1 0.1 0.1 0.1 Anderson 2007 [8] η0.2 0.2 0.2 0.2 Millard 2017 [19] fParameters Lumbar Hip Knee Ankle fA h0.210 0.720 1.042 -0.275 k1 1 1 1 p-0.173 -0.380 -0.270 -0.089 fP E kPE 3 3 3 3 εM 00.5 0.5 0.5 0.5 l00.19 0.70 1.40 -0.10 Development of MTGs for OCT of human walking Page 30 4.2 2D HAT tracking problem The main purpose of this section is to faithfully reproduce a gait cycle by using a 2D HAT model when muscle torque generators are added to it. The more similar the coordinates and the torques are to experimentally obtained data, the better the tracking and therefore, the more reliable the model. In order to quantify how well the tracking is performed, the Root Mean Square Error (RMSE) (see Equation 9) has been taken into account throughout this section. As it is stated in section 3.2.3, several studies are carried out in order to determine the final MTG version that will eventually be compared to the torque-driven model. Therefore, in this section it is discussed the process that leads to this final version as well as the comparison between the two models. 4.2.1 Study of different formulations in the MTG cost function In a first attempt to choose the cost function that yields to best results in terms of tracking, we calculated the RMSE of coordinates (see Table 7) and torques (see Table 8) at each joint for the different cost functions A, B, C and Dpreviously defined in section 3.2.3. Also note that the RMSE mean for each cost function is computed in both cases. Table 7: RMSE of coordinates obtained from the different formulations. Cost function Lumbar extension [◦] Hip flexion [◦]Knee angle [◦]Ankle angle [◦] Right Left Right Left Right Left Mean [◦] A 0.81 2.66 1.83 1.86 0.91 3.16 8.84 2.87 B 1.12 3.39 2.82 2.60 1.98 3.23 6.43 3.08 C 0.73 2.82 1.74 2.02 1.12 3.23 6.58 2.60 D 0.79 3.43 2.11 3.09 2.69 3.66 6.76 3.22 Table 8: RMSE of torques obtained from the different formulations. Cost function Lumbar extension [Nm] Hip flexion [Nm]Knee angle [Nm]Ankle angle [Nm] Right Left Right Left Right Left Mean [Nm] A 8.97 10.41 8.84 4.77 5.17 2.17 3.04 6.20 B 5.84 7.22 6.33 3.89 3.70 2.10 3.02 4.59 C 8.53 10.62 9.74 4.80 5.44 2.16 3.33 6.37 D 7.61 10.33 9.76 5.37 5.70 1.55 3.23 6.22 Clearly, the cost function that provides best results is Bas it shows significantly lower RMSE values, specially when tracking torques. The root cause of this value is that cost function Bis precisely the one that has a highest weighted factor in the term where torques are involved. That leads to a higher penalization if the difference between the computed torque and the experimental one is higher as well. Bearing in mind that the differences when tracking coordinates are not that significant between the different cost functions, it was decided to keep B(see Equation 10) as the cost function in the next study. J= nq X i=1 (qi−qexpi)2+ nq X i=1 (τi−τexpi)2+ 0.001 nq X i=1 ˙a2 i+ 0.01 nq X i=1 ¨q2 i(10) Development of MTGs for OCT of human walking Page 31 4.2.2 Influence of passive elements and joint damping in motion tracking Similarly to the previous study, we compute the RMSE to see whether adding passive elements and the joint damping to the model can indeed lead to more realistic results (see Table 9and Table 10). Remember that these elements are added to supress possible vibrations in the musculature and thus are supposed to give a better insight into the human walking behaviour. Table 9: RMSE of coordinates with and without passive elements (PE) and joint damping (DAMP). DAMP & PE Lumbar extension [◦] Hip flexion [◦]Knee angle [◦]Ankle angle [◦] Right Left Right Left Right Left Mean [◦] With 1.79 1.87 3.12 2.06 6.92 2.99 5.20 3.42 Without 1.12 3.39 2.82 2.60 1.98 3.23 6.43 3.08 Table 10: RMSE of torques with and without passive elements (PE) and joint damping (DAMP). DAMP & PE Lumbar extension [Nm] Hip flexion [Nm]Knee angle [Nm]Ankle angle [Nm] Right Left Right Left Right Left Mean [Nm] With 7.27 6.98 7.21 3.21 3.72 1.98 2.93 4.76 Without 5.84 7.22 6.33 3.89 3.70 2.10 3.02 4.59 In this case, the final choice is not that clear. While more than half of the joints present lower RMSE values when these elements are added in both cases (coordinates and torques), it is still not enough to the overall contribution. That is due to the fact that in some cases, the RMSE calculated are notably higher compared to the version without these elements. Some examples are the left knee angle in coordinates or the lumbar extension in torques. If a further exploration of the parameters defining these elements was conducted, that would probably yield lower errors and hence, a better tracking. Despite the fact that a formulation without these elements seems to provide better results, a comparison between both computational times was conducted in order to secure a final choice. The results were that the computational time with the passive elements and joint damping was more than twice as long as without them. Therefore, the final study was carried out without these elements. 4.2.3 Evaluation of muscle torque generators implementation Once the final version of the MTG model is clearly defined, we can now proceed to the final study. As it was previously stated, in this last section, the MTG model is compared to the torquedriven model, which is a previous version of the 2D HAT model at skeletal level. Likewise the previous studies, the RMSE of coordinates (see Table 11) and torques (see Table 12) at each joint will also be of interest in this final study. Table 11: RMSE of coordinates in the torque-driven model (TD) and the MTG model. Model Lumbar extension [◦] Hip flexion [◦]Knee angle [◦]Ankle angle [◦] Right Left Right Left Right Left Mean [◦] TD 0.41 0.47 0.46 0.64 0.61 0.84 0.70 0.59 MTG 1.12 3.39 2.82 2.60 1.98 3.23 6.43 3.08 Development of MTGs for OCT of human walking Page 32 Table 12: RMSE of torques in the torque-driven model (TD) and the MTG model. Model Lumbar extension [Nm] Hip flexion [Nm]Knee angle [Nm]Ankle angle [Nm] Right Left Right Left Right Left Mean [Nm] TD 3.95 4.54 6.15 2.49 2.16 1.79 1.94 3.29 MTG 5.84 7.22 6.33 3.89 3.70 2.10 3.02 4.59 When computing the RMSE in the final study, we obtain considerably low values in both models, which means that the tracking is well performed in both cases. More specifically, if we take a look at the obtained errors in torques, we note that in general the difference between models is quite low. For instance, the obtained RMSE for the left hip flexion torque in the torque-driven model is 6.15 Nm and similarly, a value of 6.33 Nm in the MTG version is obtained. Nonetheless, in some joints the difference between models is still relevant. For instance, the joint in the MTG model that presents the maximum relative error obtained when tracking coordinates is the left ankle with an RMSE value of 6.43ocompared to 0.70oin the torque-driven version. Since the RMSE only takes into account the average error between functions, it is decided to complement this information with an evolution of the coordinates and torques for each joint. If we take a look at Figure 21 a more detailed comparison between the tracking motion in both models is carried out. Therefore, it is possible to check if the tracked results show a similar tendency to the experimental data. In order to ease the visualization of the results, only left joints were plotted. We decided to show left joints instead of the right ones because the stance phase begins with the left leg. Figure 21: Tracking of coordinates and torques in a 2D HAT model based in a torque-driven (in blue) and an MTG (in red) versions. Experimental data is shown with a dashed black line. The gait accounts the 80% of the cycle and left joints are shown. Development of MTGs for OCT of human walking Page 33 If we take a look at the hip left flexion again, it is possible to see that despite the fact that RMSE in torques is similar between models, the error is still a bit high compared, for example, to the tracked torque in the ankle angle. The hip flexion angle shows much better results than its torque version. Note that when MTGs are considered in coordinates, a slight difference in the peak value is produced. Nontheless, the rest of the tracking is generally adequate. Moving now to the knee angle, similar results to the hip flexion are obtained. For instance, there is also a small difference between the peak values in the tracked coordinates when MTGs are employed, but the rest of the tracking is fine. If we now analyse the torques, the torquedriven model seems that is performing a better tracking. That can also be explained through the RMSE obtained. While the value obtained with the torque-driven version is only 2.16 Nm, when we apply MTGs a value of 3.70 Nm is reported. As it was previously stated, the ankle angle shows the maximum relative error in terms of coordinate tracking when MTGs are applied. However, the tendency of it is still adequate. Moving now to the torque tracking analysis, we observe an almost perfect tracking in both models. This notorious difference between coordinates and torques suggests that further exploration regarding the MTG implementation in the ankle should be undertaken. Finally, if we give some insight into the observed results in the lumbar extension angle, the torque-driven model provides better results again, especially at the beginning of the motion. Nonetheless, the tendency of the tracking is correct and generally well adjusted to the experimentally obtained data. Conversely, tracking in torques are a bit further from the experimental values. Similarly to the hip flexion case, both models perform a similar tracking while none of them is perfectly adjusted to the experimental data. In order to give a deeper insight into the pairs of muscle torque generators involved in each joint, we decided to break them down into their flexor and extensor part. Figure 22 shows the contribution of each MTG pair throughout the cycle. Figure 22: Pairs of muscle torque generators at each left joint and lumbar. Development of MTGs for OCT of human walking Page 34 Two remarkable facts can be seen clearly seen in Figure 22. On the one hand, the different MTG pairs are almost symmetric to one another, except in the left ankle at the beginning of the gait cycle, in which a higher contribution of the extensor torque is reported. On the other hand, we obtained quite high unexpected values in some parts of the cycle in each of the joints. Bearing in mind how these pairs of MTGs are calculated in the final version of the MTG model (see Equation 11), all terms are the same except the isometric maximum torque τoand the muscle activations a. τM i=τM oafA(q)fV( ˙q)(11) Therefore, a further exploration in these two parameters could be conducted in order to gain better insight into the contribution of each the flexor and extensor torques in the gait cycle. 4.2.4 General discussion By and large, the tracking performed when muscle torque generators are taken into account is adequate since no extreme RMSE are reported nor any undesired tendency in the tracking is observed for the plotted joints. It is also remarkable that a better tracking in the torques is achieved in comparison to the coordinates. However, the tracking is still better performed by the torque-driven model, especially when tracking coordinates. That is probably due to the fact that the cost function considered in the MTG-driven model penalizes more the torque errors rather than the coordinate ones. Therefore, if we wanted to obtain better results in terms of tracking coordinates new readjustments in the cost function should be considered. Another fact that has not yet been mentioned in this section is the computational time as well as the iterations needed to obtain the optimal solution. Table 13 shows both the time and the number of iterations that the solver needed to converge to an optimal solution in the studied tracking problem. Table 13: Computational time and number of iterations required in both models Time [s] Iterations TD 2.91 23 MTG 92.87 1403 As it was expected, the model in which muscle torques generators are employed shows higher values in both the computational time and number of iterations. That is because more restrictions have to be satisfied in the optimal control formulation regarding the torques in comparison to the torque-driven version. Finally, despite the fact that the tracking performed by the MTG model is correct, it was expected that the model had provided a more realistic behaviour in comparison to the torquedriven model. That is due to the fact that still further explorations need to be conducted in some of the parameters that define the pairs of MTGs before this model is employed for other purposes such human motion prediction. Parameters such as the maximum isometric torque or the characteristic musculotendon curves involved in each pair of the MTGs could be more in depth investigated to gain better insights. Development of MTGs for OCT of human walking Page 35 5 Project impact In this section we provide a detailed explanation about the economic cost of the project. Besides, we contemplate the environmental impact produced in this work. . 5.1 Economic cost The economic cost of the project is based on four different aspects: the depreciation of the computer used for the simulations, the MATLAB and GPOPS-II licenses, the working time (student and supervisors), and finally, the electrical energy consumed. Depreciation is calculated from the total price of the electronic device, its useful time and the total time it has been used. The useful life of the computer is considered to be 4 years, which corresponds to a total of 35040 hours. Finally, the laptop has been used for an amount of 320 hours during the thesis. MATLAB and GPOPS-II licenses expire after one year. Therefore, their useful life is about 8760 hours. Considering that the software was used only during half of the hours when the computer was used, it is then estimated that both MATLAB and GPOPS-II have been used for a total of 160 hours. Regarding the student’s hours of dedication, it should be taken into account the time in which the computer was used (320 hours), plus 30 hours of meetings, calculations and reflections. The hours of support and supervision are considered to be 40 hours in total (considering both supervisors). Finally, the electrical power of light and computers have been estimated to be 35 W for each of them. A light bulb was only switched on during the hours of less light (a total of 100 hours) and the computer during the previously mentioned 320 hours. We assumed constant the price of electricity, with a value of 0.29 €/kWh. As can be seen in Table 14, the total cost of the project is 3428.28 €. Table 14: Calculation of the final project cost. Variable costs of the computer and licenses are obtained from dividing their fixed cost by their useful life in hours. Variable cost of electrical energy is found from its price multiplied by the power consumption of computers and light. Fixed cost [€]Useful life [years] Variable cost [€/h] Time spent in the project [h] Cost related to the ptoject [€] Laptop 750 4 0.0214 320 6.85 MATLAB license 840 1 0.0959 160 15.34 GPOPS-II license 100 1 0.0114 160 1.83 Student - - 8 350 2800.00 Supervisors - - 15 40 600.00 Electrical energy - - 0.0102 420 4.26 TOTAL COST 3428.28 € Development of MTGs for OCT of human walking Page 36 5.2 Environmental impact The environmental impact of this project is minimal since the main studies were based on simulations using a computer at home. However, electricity and electronic devices such as a laptop and a tablet can be considered. Regarding the electrical consumption, a computer and a light bulb have been assumed to be switched on while working on the thesis. However, their consumption is marginal and thus, they are not considered to have such a negative environmental impact in the project. When evaluating the deterioration of the electronic devices we have to take into account that once their useful life is finished, they must be treated individually. This evaluation is in agreement with the Regulation 2017/699 [34], which establishes a common methodology for the calculation of the weight of electrical and electronic equipment (EEE) as well as for the calculation of the quantity of waste in electrical and electronic equipment (WEE). Development of MTGs for OCT of human walking Page 37 Conclusions In this bachelor’s thesis, optimal control techniques have been applied to both the prediction of a new motion and the tracking of a known one. In order to implement the optimal control algorithm, we used GPOP S −II, which works in MATLAB. Therefore, there has been a process of getting acquainted with optimal control theory and its environment. In both formulations, a torque-driven model (i.e., at skeletal level) is obtained from the OpenSim repository. Subsequently, we employed muscle torque generators with MATLAB in each of them as an alternative of detailed musculoskeletal models. Regarding the predictive problem when using a simple pendulum model, muscle torque generators were slightly advantageous than its torque-driven version when the studied simulation encompassed a shorter range of motion. Conversely, when using a 2D HAT model for human gait tracking, despite the fact that the tracking performed with MTGs was adequate, a more realistic behaviour than its torque-driven version was expected. This project can be considered a first research work in muscle torque generators modelling in the Biomechanical Engineering Lab (BIOMEC) at UPC. Thus, further studies would need to be conducted in order to obtain more realistic models. Parameters such as the maximum isometric torque or the characteristic musculotendon curves involved in each pair of the MTGs could be more in depth investigated to gain better insights. Development of MTGs for OCT of human walking Page 38 Acknowledgements I would like to express my sincere thanks to my supervisors, Míriam Febrer Nafría and Josep María Font Llagunes, for their support, advice and for giving me the opportunity to work for the first time in a research group such as the BIOMEC Lab. I want to give thanks especially to you, Míriam, for the weekly meetings, for your constant feedback until the end of the project, and for being able to highlight the well done work as well. Moreover, I would like to thank family and friends for supporting me when things did not come that easily and for always showing interest in each of the small (or big) advances in the project. Development of MTGs for OCT of human walking Page 39 References [1] Marko Ackermann and Werner Schiehlen. “Dynamic Analysis of Human Gait Disorder and Metabolical Cost Estimation”. In: Archive of Applied Mechanics 75 (Sept. 2006), pp. 569–594. doi:10.1007/s00419-006-0027-7. [2] A. Alamdari and V.N. Krovi. “Chapter Two - A Review of Computational Musculoskeletal Analysis of Human Lower Extremities”. In: Human Modelling for Bio-Inspired Robotics. Ed. by Jun Ueda and Yuichi Kurita. Academic Press, 2017, pp. 37–73. isbn: 978-0-12-8031377. doi:https://doi.org / 10.1016/B978012 - 8031377.000033.url:https: //www.sciencedirect.com/science/article/pii/B9780128031377000033. [3] BBC. Muscular system. Agonist and antagonist muscle pairs. url:https://www.bbc.co.uk/ bitesize/guides/z8stfrd/revision/4. (accessed: 15.06.2022). [4] Peter Brown and John McPhee. “A 3D ellipsoidal volumetric foot–ground contact model for forward dynamics”. In: Multibody System Dynamics 42 (Apr. 2018). doi:10.1007 / s11044-017-9605-4. [5] Jeremy C O’Connor Christopher L Vaughan Brian L Davis. Dynamics of Human Gait. Kiboho Publishers, 1992. isbn: 0-620-23558-6. [6] David Civantos. “Study of the foot-ground contact model for the prediction of human gait”. In: (July 2020). doi:https://upcommons.upc.edu/handle/2117/349395. [7] Matthew Millard Thomas Uchida Ajay Seth Scott L. Delp. “Flexing Computational Muscle: Modeling and Simulation of Musculotendon Dynamics”. In: Journal of Biomechanical Engineering 135 (Feb. 2013), pp. 70–77. doi:https://doi.org/10.1115/1.4023390. [8] Michael L. Madigan Dennis E. Anderson and Maury A. Nussbaum. “Maximum voluntary joint torque as a function of joint angle and angular velocity: Model development and application to the lower limb”. In: Journal of Biomechanics 40 (Mar. 2007), pp. 3105– 3113. url:https://www.sciencedirect.com/journal/journal-of-biomechanics. [9] OpenSim Documentation. Characteristic Musculotendon Curves.url:https : / / simtk - confluence.stanford.edu:8443/display/OpenSim/Characteristic+Musculotendon+ Curves. (accessed: 31.05.2022). [10] Marcus G. Pandy Frank C. Anderson. “Dynamic Optimization of Human Walking”. In: Journal of Biomechanical Engineering 123 (Oct. 2001), pp. 381–390. doi:10.1115/1.1392310. url:chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://homes.cs. washington.edu/~todorov/courses/amath533/AndersonPandy.pdf. [11] Benjamin Fregly. “Design of Optimal Treatments for Neuromusculoskeletal Disorders using Patient-Specific Multibody Dynamic Models”. In: International journal for computational vision and biomechnanics 2 (July 2009), pp. 145–155. [12] Benjamin J. Fregly. “A Conceptual Blueprint for Making Neuromusculoskeletal Models Clinically Useful”. In: Applied Sciences 11.5 (2021). issn: 2076-3417. doi:10 . 3390 / app11052037.url:https://www.mdpi.com/2076-3417/11/5/2037. [13] Friedl De Groote - Allison L Kinney - Anil V Rao - Benjamin J Fregly. “Evaluation of Direct Collocation Optimal Control Problem Formulations for Solving the Muscle Redundancy Problem”. In: Epub 44(10) (Oct. 2016), pp. 2922–2936. doi:10.1007/s10439-016-1591-9. [14] Archibald Vivian Hill. “The heat of shortening and the dynamic constants of muscle”. In: Proc. R. Soc. Lond. 126 (Oct. 1938). doi:B126136âĂŞ195. [15] McNally W Inkol K.A. Brown C. “Muscle torque generators in multibody dynamic simulations of optimal sports performance”. In: Multibody System Dynamics 50 (Dec. 2020), pp. 435–452. doi:10.1007/s11044-020-09747-9.url:https://link.springer.com/ article/10.1007/s11044-020-09747-9.