Manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 13
Manual SimTK optcntrlmuscle (v2.1)
Friedl De Groote, Maarten Afschrift, Tom Van Wouwe, Antoine Falisse
07/11/2018
Contents
1 Release notes 2
1.1 Release2.1 ....................................... 2
2 Overview 2
3 Installation Instruction 3
4 Main Function 3
4.1 UsingGPOPS...................................... 3
4.1.1 With explicit activation dynamics formulation (De Groote et al. (2016)) . 3
4.1.2 With implicit activation dynamics formulation (De Groote et al. (2009)) . 3
4.2 UsingCasADi...................................... 4
4.2.1 With explicit activation dynamics formulation (De Groote et al. (2016)) . 4
4.2.2 With implicit activation dynamics formulation (De Groote et al. (2009)) . 4
4.3 InputArguments.................................... 4
4.4 Outputarguments ................................... 5
4.4.1 UsingGPOPS ................................. 5
4.4.2 UsingCasADi.................................. 6
5 GPOPS-II 7
5.1 Setup .......................................... 7
5.2 Output ......................................... 7
6 CasADi 8
6.1 Setup .......................................... 8
6.2 Output ......................................... 8
7 Muscle model 8
8 Tips and Tricks 8
9 Examples 8
9.1 Walking example De Groote et al. 2016 . . . . . . . . . . . . . . . . . . . . . . . 8
9.2 Running example De Groote et al. 2016 . . . . . . . . . . . . . . . . . . . . . . . 10
9.3 OpenSim installation example Gait10dof18m . . . . . . . . . . . . . . . . . . . . 11
9.4 OpenSim installation example Gait23dof54m . . . . . . . . . . . . . . . . . . . . 12
1
1 Release notes
1.1 Release 2.1
CasADi was added as an alternative to GPOPS-II and ADiGator (see section 2 for details).
The reserve actuators (RActivation) were unscaled in the output of the main functions.
The time derivatives of the muscle contraction dynamics states, i.e. normalized muscle
fiber velocities or derivatives of normalized tendon forces, were added to the cost function
with a small weighting factor to prevent spiky outputs.
The tendon stiffness was added as an optional user parameter (see section 9.3 for example).
2 Overview
The provided MATLAB code solves the muscle redundancy problem using direct collocation as
described in De Groote F, Kinney AL, Rao AV, Fregly BJ. Evaluation of direct collocation opti-
mal control problem formulations for solving the muscle redundancy problem. Annals of Biomedi-
cal Engineering (2016). http://link.springer.com/article/10.1007%2Fs10439-016-1591-9.
From v2.1, CasADi can be used as an alternative to GPOPS-II and ADiGator. CasADi is
an open-source tool for nonlinear optimization and algorithmic differentiation (https://web.
casadi.org/). Results using CasADi and GPOPS-II are very similar (differences can be at-
tributed to the different direct collocation formulations and scaling). We used CasADi’s Opti
stack, which is a collection of CasADi helper classes that provides a close correspondence
between mathematical NLP notation and computer code (https://web.casadi.org/docs/
#document-opti). CasADi is actively maintained and developed, and has an active forum
(https://groups.google.com/forum/#!forum/casadi-users).
From v1.1, an implicit formulation of activation dynamics can be used to solve the muscle re-
dundancy problem. Additionally, by using the activation dynamics model proposed by Raasch
et al. (1997), we could introduce a nonlinear change of variables to exactly impose activa-
tion dynamics in a continuously differentiable form, omitting the need for a smooth approx-
imation such as described in De Groote et al. (2016). A result of this change of variables
is that muscle excitations are not directly accessible during the optimization. Therefore, we
replaced muscle excitations by muscle activations in the objective function. This implicit for-
mulation is described in De Groote F, Pipeleers G, Jonkers I, Demeulenaere B, Patten C,
Swevers J, De Schutter J. A physiology based inverse dynamic analysis of human gait: potential
and perspectives F. Computer Methods in Biomechanics and Biomedical Engineering (2009).
http://www.tandfonline.com/doi/full/10.1080/10255840902788587. Results from both
formulations are very similar (differences can be attributed to the slightly different activation
dynamics models and cost functions). However, the formulation with implicit activation dy-
namics (De Groote et al., (2009)) is computationally faster. This can mainly be explained by
the omission of a tanh function in the constraint definition, whose evaluation is computationally
expensive when solving the NLP.
Any questions? Please contact us: friedl.degroote@kuleuven.be,maarten.afschrift@kuleuven.be,
tom.vanwouwe@kuleuven.be, and antoine.falisse@kuleuven.be.
2
3 Installation Instruction
Add the main folder and subfolder to your MATLAB path
1addpath(genpath('C/......./SimTK optcntrlmuscle'))).
Several software packages are needed to run the program
The OpenSim MATLAB interface is used to generate the inputs to the optimal control
problem based on a scaled OpenSim model and the solution of inverse kinematics (pro-
viding the solution of inverse dynamics is optional). To this aim, install OpenSim and set
up the OpenSim MATLAB interface (OpenSim: https://simtk.org/frs/?group_id=
91, OpenSim API: http://simtk-confluence.stanford.edu:8080/display/OpenSim/
Scripting+with+Matlab).
Using GPOPS
–GPOPS-II is used to solve the optimal control problem using direct collocation (http:
//www.gpops2.com/). A one-time 30-day trial license is avaiable for all users who
register.
–ADiGator is used for automatic differentiation (https://sourceforge.net/projects/
adigator/).
Using CasADi
–CasADi is used for nonlinear optimization and algorithmic differentiation (https:
//web.casadi.org/).
4 Main Function
SolveMuscleRedundancy is the main function of this program and is used to solve the muscle
redundancy problem. There are eight variants of this function that differ based on the chosen
optimization package (GPOPS or CasADi), the formulation of the contraction dynamics (nor-
malized tendon force or normalized muscle fiber length as a state), and the formulation of the
activation dynamics (explicit or implicit).
4.1 Using GPOPS
4.1.1 With explicit activation dynamics formulation (De Groote et al. (2016))
SolveMuscleRedundancy FtildeState GPOPS uses the normalized tendon force as a state
SolveMuscleRedundancy lMtildeState GPOPS uses the normalized muscle fiber length as
a state
4.1.2 With implicit activation dynamics formulation (De Groote et al. (2009))
SolveMuscleRedundancy FtildeState actdyn GPOPS uses the normalized tendon force as
a state
SolveMuscleRedundancy lMtildeState actdyn GPOPS uses the normalized muscle fiber
length as a state
3
4.2 Using CasADi
4.2.1 With explicit activation dynamics formulation (De Groote et al. (2016))
SolveMuscleRedundancy FtildeState CasADi uses the normalized tendon force as a state
SolveMuscleRedundancy lMtildeState CasADi uses the normalized muscle fiber length as
a state
4.2.2 With implicit activation dynamics formulation (De Groote et al. (2009))
SolveMuscleRedundancy FtildeState actdyn CasADi uses the normalized tendon force as
a state
SolveMuscleRedundancy lMtildeState actdyn CasADi uses the normalized muscle fiber
length as a state
4.3 Input Arguments
Required input arguments for SolveMuscleRedundancy
1. model path: directory and filename of the scaled OpenSim model (.osim file). The code
should work with any OpenSim model with valid muscle-tendon parameters for which
OpenSim’s Inverse Dynamics and Muscle Analysis Tools generate reliable results. Note
that only the muscle-tendon parameters and not the muscle model specified in the osim-file
are used (for details see Muscle model).
2. IK path: directory and filename of the inverse kinematics solution (.mot file).
3. ID path: directory and filename of the inverse dynamics solution (.sto file). If left empty,
the inverse dynamics solution will be computed from the external loads (see Optional input
arguments).
4. time: 1 x 2 MATLAB array with the initial and final time of the analysis in seconds.
Initial and final states influence the optimal controls over a period of about 50 ms at the
beginning and end of the time interval over which the optimal control problem is solved.
Since in practice the initial and final states are generally unknown, problems should be
solved for a time interval containing five additional data points (considering a 100Hz
sampling frequency) at the beginning and end of the motion cycle. Those additional data
points should not be considered in further analyses. The user should thus not be surprised
to observe unrealistically high muscle activation at the beginning of the motion (more
details in companion paper).
5. Out path: directory where you want to store the results from the muscle analysis.
6. Misc: miscellaneous input arguments
DofNames Input is a cell array specifying for which degrees of freedom you want to
solve the muscle redundancy problem. Typically the muscle redundancy problem is
solved for one leg at a time (there are no muscles spanning both legs).
MuscleNames Input is a cell array that specifies the muscles to be included when
solving the muscle redundancy problem. All muscles that actuate (i.e. have a moment
arm with respect to) the degrees of freedom specified in DofNames Input will be
selected by default if this array is left empty.
4
Optional input arguments for SolveMuscleRedundancy
1. Misc.Loads path: directory and filename of the external loads (.xml file). The program
will use the OpenSim libraries to solve the inverse dynamics problem when the required
input argument ID path is empty and Misc.Loads path points to an external loads file.
2. Misc.ID ResultsPath: directory where the inverse dynamics results will be saved when
the input argument ID path is left empty.
3. Misc.f cutoff ID: cutoff frequency for the butterworth recursive low pass filter applied to
the inverse dynamics data (default is 6 Hz).
4. Misc.f order ID: order of the butterworth recursive low pass filter applied to the inverse
dynamics data (default is 6).
5. Misc.f cutoff LMT : cutoff frequency for the butterworth recursive low pass filter applied
to the muscle tendon lengths from the muscle analysis (default is 6 Hz).
6. Misc.f order LMT : order of the butterworth recursive low pass filter applied to the muscle
tendon lengths from the muscle analysis (default is 6).
7. Misc.f cutoff dM : cutoff frequency for the butterworth recursive low pass filter applied to
the muscle moment arms from the muscle analysis (default is 6 Hz).
8. Misc.f order dM : order of the butterworth recursive low pass filter applied to the muscle
moment arms from the muscle analysis (default is 6).
9. Misc.f cutoff IK : cutoff frequency for the butterworth recursive low pass filter applied
to the inverse kinematics data (default is 6 Hz) when performing the muscle analysis to
compute muscle-tendon lengths and moment arms.
10. Misc.f order IK : order of the butterworth recursive low pass filter applied to the inverse
kinematics data (default is 6).
11. Misc.Mesh Frequency: number of mesh interval per second (default is 100, but a denser
mesh might be required to obtain the desired accuracy especially for faster motions).
12. Misc.Atendon: vector with tendon stiffness for the selected muscles. The order should
correspond to MuscleNames Input. The default value is 35 and a lower value corresponds
to a more compliant tendon. The default value will be used when left empty. An example
is provided in section 9.3 to set a different stiffness to the Achilles tendon.
4.4 Output arguments
4.4.1 Using GPOPS
1. Time: time vector.
2. MExcitation: optimal muscle excitation (matrix dimension: number of collocation points
x number of muscles).
3. MActivation: optimal muscle activation (matrix dimension: number of collocation points
x number of muscles).
5
4. RActivation: activation of the reserve actuators (matrix dimension: number of collocation
points x number of degrees of freedom).
5. TForcetilde: normalized tendon force (matrix dimension: number of collocation points x
number of muscles).
6. TForce: tendon force (matrix dimension: number of collocation points x number of mus-
cles).
7. lMtilde: normalized muscle fiber length (matrix dimension: number of collocation points
x number of muscles).
8. lM: muscle fiber length (matrix dimension: number of collocation points x number of
muscles) .
9. MuscleNames: cell array that contains the names of the selected muscles (matrix dimen-
sion: number of muscles).
10. OptInfo: output structure created by GPOPS-II.
11. DatStore: data structure with input information for the optimal control problem.
4.4.2 Using CasADi
CasADi uses piecewise-constant controls in the mesh intervals. We therefore distinguish between
collocation points (for the states) and mesh points (for the states and controls) in the output
arguments.
1. Time: time vector
(a) Time.meshPoints: time at mesh points
(b) Time.collocationPoints: time at collocation points
2. MExcitation.meshPoints: muscle excitation (matrix dimension: number of mesh points
points x number of muscles).
3. MActivation: muscle activation
(a) MActivation.meshPoints: muscle activation at mesh points (matrix dimension: num-
ber of mesh points x number of muscles).
(b) MActivation.collocationPoints: muscle activation at collocation points (matrix di-
mension: number of collocation points x number of muscles).
4. RActivation.meshPoints: activation of the reserve actuators (matrix dimension: number
of mesh points x number of degrees of freedom).
5. TForcetilde: normalized tendon force
(a) TForcetilde.meshPoints: normalized tendon force at mesh points (matrix dimension:
number of mesh points x number of muscles).
(b) TForcetilde.collocationPoints: normalized tendon force at collocation points (matrix
dimension: number of collocation points x number of muscles).
6. TForce: tendon force
6
(a) TForce.meshPoints: tendon force at mesh points (matrix dimension: number of mesh
points x number of muscles).
(b) TForce.collocationPoints: tendon force at collocation points (matrix dimension: num-
ber of collocation points x number of muscles).
7. lMtilde: normalized muscle fiber length
(a) lMtilde.meshPoints: normalized muscle fiber length at mesh points (matrix dimen-
sion: number of mesh points x number of muscles).
(b) lMtilde.collocationPoints: normalized muscle fiber length at collocation points (ma-
trix dimension: number of collocation points x number of muscles).
8. lM: muscle fiber length
(a) lM.meshPoints: muscle fiber length at mesh points (matrix dimension: number of
mesh points x number of muscles).
(b) lM.collocationPoints: muscle fiber length at collocation points (matrix dimension:
number of collocation points x number of muscles).
9. MuscleNames: cell array that contains the names of the selected muscles (matrix dimen-
sion: number of muscles).
10. OptInfo: output structure with settings used in CasADi.
11. DatStore: data structure with input information for the optimal control problem.
5 GPOPS-II
5.1 Setup
The GPOPS-II setup is accessible through the function SolveMuscleRedundancy < ... > -
GPOPS.m under GPOPS setup. The user is referred to the GPOPS-II user guide for setup
options. A higher accuracy can be reached by adjusting, for instance, the number of mesh in-
tervals. This however comes at the expense of the computational time. 100 mesh intervals per
second are used by default.
5.2 Output
The output variable OptInfo contains all information related to the optimal control problem
solution. Convergence to an optimal solution is reached when output.result.nlpinfo is flagged 0
(”EXIT: Optimal solution found” in the command window of MATLAB). The mesh accuracy
can be assessed with output.result.maxerrors. Cost functional, control, state (and costate) can
be accessed in output.result.solution.phase.
To recall, the user should consider extending the time interval by 50-100 ms at the beginning and
end of the motion to limit the influence of the unknown initial and final state on the solution.
Results from those additional periods should not be considered realistic and will typically result
in high muscle activation.
7
6 CasADi
6.1 Setup
The CasADi setup is accessible through the function SolveMuscleRedundancy < ... > CasADi.m
under CasADi setup. The user is referred to the CasADi user guide for setup options.
6.2 Output
The output variable OptInfo contains information related to the optimal control problem set-
tings. As with GPOPS, the user should consider extending the time interval by 50-100 ms at
the beginning and end of the motion to limit the influence of the unknown initial and final state
on the solution. Results from those additional periods should not be considered realistic and
will typically result in high muscle activation.
7 Muscle model
The musculotendon properties are fully described in the supplementary materials of the afore-
mentioned publication. Importantly, only the tendon slack length, optimal muscle fiber length,
maximal isometric muscle force, optimal pennation angle and maximal muscle fiber contraction
velocity are extracted from the referred OpenSim model. Other properties are defined in the code
and can be changed if desired. By default, the activation and deactivation time constants are
15 and 60 ms respectively (see tau act and tau deact in SolveMuscleRedundancy < state >.m).
8 Tips and Tricks
1. We advise users to use the formulations with normalized tendon force as a state. They
appear to have improved convergence properties as compared to formulations with nor-
malized muscle fiber length as a state. This might be due to the more linear relationship
between muscle activation and tendon force than between muscle activation and muscle
length.
2. We advise users to use the formulations with implicit activation dynamics as they are in
most cases computationally more efficient.
3. If you observe spiky output variables (e.g. muscle excitations), try increasing the mesh fre-
quency. This might improve the results although it might also increase the computational
time.
9 Examples
Four examples are provided in the folder examples.
9.1 Walking example De Groote et al. 2016
1 c l e a r a l l ; c l o s e a l l ; c l c
2
3%% Choose o p t i m i z a t i o n f ra mewor k
4% fra m ework = GPOPS ;
5 framework = CasADi ;
6
8
7%% Choose c o n t r a c t i o n dynami cs f o r mu l a t io n
8% f o r m u l a t i o n c o n t d y n = l M t i l d e S t a t e ;
9 f o r m u l a t i o n c o n t d y n = FtildeState ;
10
11 %% Choose a c t i v a t i o n dynamics f o r m u l a t io n
12 % f o r m u l a t i o n a c t d y n = DeG root e20 16 ;
13 formulation actdyn = DeGroote2009 ;
14
15 %% Example
16 % Add main f o l d e r and s u b f o l d e r t o matla b path ( i n s t a l l a t i o n )
17 f i l e p a t h=which ( Walking DeGrooteetal2016.m ) ;
18 [ DirExample Walking ,¬,¬]= f i l e p a r t s ( f i l e p a t h ) ;
19 [ DirExample ,¬]= f i l e p a r t s ( D i rExample Walking ) ;
20 [ MainD ir , ¬]= f i l e p a r t s ( DirExample ) ;
21 add path ( ge n p a t h ( MainDir ) ) ;
22
23 % Needed I np u t Arguments
24 I K p a th= f u l l f i l e ( MainDir , Examples ,Walking DeGrooteetal2016 ,WalkingData ,Walking IK.mot ) ;
25 I D p ath= f u l l f i l e ( MainDir , Examples ,Walking DeGrooteetal2016 ,WalkingData ,Walking ID.sto ) ;
26 m od e l pa t h= f u l l f i l e ( MainDir , Examples ,Walking DeGrooteetal2016 ,WalkingData ,subject1.osim ) ;
27 t im e =[0 . 5 1 6 1 . 9 5 ] ; % R i gh t s t a n c e p has e (+50ms b e g i n n i n g and e nd o f t ime i n t e r v a l , more ...
d e t a i l s s e e manual and p u b l i c a t i o n )
28 Out path= f u l l f i l e ( MainDir , Examples ,Walking DeGrooteetal2016 ,Results ) ;
29
30 Misc.DofNames Input={ankle angle r ,k n e e a n g l e r ,h i p f l e x i o n r ,hip rotation r ,h i p a d d u c t i o n r };
31 Misc.MuscleNames Input={};% S e l e c t s a l l m usc le s f o r i np u t DOFs when empty
32
33 % O p t i o n a l I n p u t Ar gu men ts
34 M i s c . f c u t o f f I D = 6 ; % c u t o f f f r e q u e n c y f i l t e r i n g ID
35 Misc.f order ID = 4; % o r d e r f r e q u e n c y f i l t e r i n g ID
36 M i s c . f c u t o f f l M T = 6 ; % c u t o f f f r e q u e n c y f i l t e r i n g lMT
37 Misc.f order lMT = 4; % o r d e r f r e q u e n c y f i l t e r i n g lMT
38 M i s c . f c u t o f f d M= 6 ; % c u t o f f f r e q u e n c y f i l t e r i n g MA
39 M i s c . f o r d e r d M = 4 ; % o r d e r f r e q u e n c y f i l t e r i n g MA
40 M i s c . f c u t o f f I K= 6 ; % c u t o f f f r e q u e n c y f i l t e r i n g IK
41 M i s c . f o r d e r I K = 4 ; % o r d e r f r e q u e n c y f i l t e r i n g IK
42
43 %% S o l v e t h e p r ob le m
44 switch framework
45 c a s e GPOPS
46 switch formulation actdyn
47 c a s e DeGroote2016
48 switch formulation contdyn
49 c a s e lMtildeState
50 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
51 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
52 SolveMuscleRedundancy lMtildeState GPOPS ( ...
53 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
54 c a s e FtildeState
55 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
56 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
57 SolveMuscleRedundancy FtildeState GPOPS ( . . .
58 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
59 end
60
61 c a s e DeGroote2009
62 switch formulation contdyn
63 c a s e lMtildeState
64 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
65 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
66 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
67 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
68 SolveMuscleRedundancy lMtildeState actdyn GPOPS ( ...
69 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
70 c a s e FtildeState
71 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
72 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
73 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
74 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
75 SolveMuscleRedundancy FtildeState actdyn GPOPS ( ...
76 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
77 end
78 end
79 c a s e CasADi
80 switch formulation actdyn
81 c a s e DeGroote2016
82 switch formulation contdyn
83 c a s e lMtildeState
84 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
85 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
86 S o l v e M u s c l e R e d u n d a n c y l M t i l d e S t a t e C a s A D i ( ...
87 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
88 c a s e FtildeState
89 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
90 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
91 SolveMuscleRedundancy FtildeState CasADi ( ...
92 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
93 end
94
95 c a s e DeGroote2009
96 switch formulation contdyn
97 c a s e lMtildeState
98 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
99 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
9
100 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
101 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
102 SolveMuscleRedundancy lMtildeState actdyn CasADi ( ...
103 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
104 case FtildeState
105 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
106 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
107 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
108 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
109 SolveMuscleRedundancy FtildeState actdyn CasADi ( . . .
110 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
111 end
112 end
113 end
9.2 Running example De Groote et al. 2016
1 c l e a r a l l ; c l o s e a l l ; c l c
2
3%% Choose o p t i m i z a t i o n f ra mewor k
4% fra m ework = GPOPS ;
5 framework = CasADi ;
6
7%% Choose c o n t r a c t i o n dynami cs f o r mu l a t io n
8% f o r m u l a t i o n c o n t d y n = l M t i l d e S t a t e ;
9 f o r m u l a t i o n c o n t d y n = FtildeState ;
10
11 %% Choose a c t i v a t i o n dynamics f o r m u l a t io n
12 % f o r m u l a t i o n a c t d y n = D eGr oote 201 6 ;
13 formulation actdyn = DeGroote2009 ;
14
15 %% Example
16 % Add main f o l d e r and s u b f o l d e r t o matla b path ( i n s t a l l a t i o n )
17 f i l e p a t h=which ( Running DeGrooteetal2016.m ) ;
18 [ D irEx am pl e R un ni ng , ¬,¬]= f i l e p a r t s ( f i l e p a t h ) ;
19 [ DirExample ,¬]= f i l e p a r t s ( DirExa m ple Runni n g ) ;
20 [ MainD ir , ¬]= f i l e p a r t s ( DirExample ) ;
21 add path ( ge n p a t h ( MainDir ) ) ;
22
23 % Needed I np u t Arguments
24 I K p ath= f u l l f i l e ( MainDir , Examples ,Running DeGrooteetal2016 ,RunningData ,Running IK.mot ) ;
25 I D p ath= f u l l f i l e ( MainDir , Examples ,Running DeGrooteetal2016 ,RunningData ,R u n n i n g I D . s t o ) ;
26 m od e l pa t h= f u l l f i l e ( MainDir , Examples ,Running DeGrooteetal2016 ,RunningData ,subject1.osim ) ;
27 t ime =[0 . 0 5 0 . 9 8 ] ; % Ri g ht s t a n c e p ha se (+50ms b eg in n in g and end o f t im e i n t e r v a l , more d e t a i l s ...
s e e manual and p u b l i c a t i o n )
28 Out p ath= f u l l f i l e ( MainDir , Examples ,Running DeGrooteetal2016 ,Results ) ;
29
30 Misc.DofNames Input={ankle angle r ,k n e e a n g l e r ,h i p f l e x i o n r ,h i p a d d u c t i o n r ,hip rotation r };
31 Misc.MuscleNames Input={};% S e l e c t s a l l m usc le s f o r i np u t DOFs when empty
32
33 % O p t i o n a l I n p u t Ar gu men ts
34 M i s c . f c u t o f f I D = 1 0 ; % c u t o f f f r e q u e n c y f i l t e r i n g ID
35 Misc.f order ID = 5; % o r d e r f r e q u e n c y f i l t e r i n g ID
36 M i s c . f c u t o f f l M T = 1 0 ; % c u t o f f f r e q u e n c y f i l t e r i n g lMT
37 Misc.f order lMT = 5; % o r d e r f r e q u e n c y f i l t e r i n g lMT
38 M i s c . f c u t o f f d M= 1 0 ; % c u t o f f f r e q u e n c y f i l t e r i n g MA
39 M i s c . f o r d e r d M = 5 ; % o r d e r f r e q u e n c y f i l t e r i n g MA
40 M i s c . f c u t o f f I K= 1 0 ; % c u t o f f f r e q u e n c y f i l t e r i n g IK
41 M i s c . f o r d e r I K = 5 ; % o r d e r f r e q u e n c y f i l t e r i n g IK
42
43 %% S o l v e t h e p r ob le m
44 switch framework
45 c a s e GPOPS
46 switch formulation actdyn
47 c a s e DeGroote2016
48 switch formulation contdyn
49 c a s e lMtildeState
50 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
51 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
52 SolveMuscleRedundancy lMtildeState GPOPS ( ...
53 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
54 c a s e FtildeState
55 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
56 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
57 SolveMuscleRedundancy FtildeState GPOPS ( . . .
58 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
59 end
60
61 c a s e DeGroote2009
62 switch formulation contdyn
63 c a s e lMtildeState
64 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
65 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
66 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
67 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
68 SolveMuscleRedundancy lMtildeState actdyn GPOPS ( ...
69 m o d e l path , IK p ath , ID p ath , tim e , Out path , Mis c ) ;
70 c a s e FtildeState
10
71 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
72 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
73 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
74 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
75 SolveMuscleRedundancy FtildeState actdyn GPOPS ( ...
76 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
77 end
78 end
79 c a s e CasADi
80 switch formulation actdyn
81 c a s e DeGroote2016
82 switch formulation contdyn
83 c a s e lMtildeState
84 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
85 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
86 S o l v e M u s c l e R e d u n d a n c y l M t i l d e S t a t e C a s A D i ( ...
87 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
88 c a s e FtildeState
89 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
90 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
91 SolveMuscleRedundancy FtildeState CasADi ( ...
92 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
93 end
94
95 c a s e DeGroote2009
96 switch formulation contdyn
97 c a s e lMtildeState
98 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
99 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
100 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
101 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
102 SolveMuscleRedundancy lMtildeState actdyn CasADi ( ...
103 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
104 case FtildeState
105 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
106 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
107 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
108 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
109 SolveMuscleRedundancy FtildeState actdyn CasADi ( . . .
110 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
111 end
112 end
113 end
9.3 OpenSim installation example Gait10dof18m
1 c l e a r a l l ; c l o s e a l l ; c l c
2
3%% Choose o p t i m i z a t i o n f ra mewor k
4% fra m ework = GPOPS ;
5 framework = CasADi ;
6
7%% Choose c o n t r a c t i o n dynami cs f o r mu l a t io n
8% f o r m u l a t i o n c o n t d y n = l M t i l d e S t a t e ;
9 f o r m u l a t i o n c o n t d y n = FtildeState ;
10
11 %% Choose a c t i v a t i o n dynamics f o r m u l a t io n
12 % f o r m u l a t i o n a c t d y n = D eGr oote 201 6 ;
13 formulation actdyn = DeGroote2009 ;
14
15 %% Example
16 % Add main f o l d e r and s u b f o l d e r t o matla b path ( i n s t a l l a t i o n )
17 f i l e p a t h=which ( Example Gait10dof18m.m ) ;
18 [ DirExample ,¬,¬]= f i l e p a r t s ( f i l e p a t h ) ;
19 [ Di rEx amp le2 , ¬,¬]= f i l e p a r t s ( DirExample ) ;
20 [ MainD ir , ¬]= f i l e p a r t s ( DirExampl e2 ) ;
21 add path ( ge n p a t h ( MainDir ) ) ;
22
23 % Needed I np u t Arguments
24 Datapath= C : \OpenSim 3 . 3 \Models\Gait10dof18musc\OutputReference ;
25 I K p a th= f u l l f i l e ( Datapath , IK ,subject01 walk IK.mot ) ;
26 I D p at h = [ ] ; % Compute ID from t he e x t e r n a l l oa d s
27 m od e l pa t h= f u l l f i l e ( Datapath , subject01.osim ) ;
28 t im e =[0 . 7 1 . 4 ] ; % P ar t o f t h e r i g h t s t a n c e p ha se
29 Out p ath= f u l l f i l e ( MainDir , Examples ,OpenSimInstallation Gait10dof18m ,Results ) ;
30
31 Misc.DofNames Input={ankle angle r ,k n e e a n g l e r ,h i p f l e x i o n r };
32 M i s c . L o a d s p a t h= f u l l f i l e ( D ata pat h , E xp e r im e n ta l D ata ,subject01 walk grf.xml ) ;
33 M i s c . I D R e s u l t s P a t h= f u l l f i l e ( Da tap at h , ID ,inversedynamics.sto ) ;
34
35 % O p t i o n a l I n p u t A rgu men ts
36 % Here i s an exa mple o f how t o a d j u s t t h e A c h i l l e s t e n d o n s t i f f n e s s .
37 % We f i r s t add t h e i n p u t a rg um en t M u sc l eN a me s I np u t w i th ALL m u s c l e s
38 % t h a t a c t u a t e t he d e g r e e s o f f re ed om l i s t e d i n Do fN ame s I np ut.
39 Misc.MuscleNames Input={hamstrings r ,b i f e m s h r ,g l u t m a x r ,. . .
40 iliopsoas r ,r e c t f e m r ,v a s t i r ,gastroc r ,soleus r ,. . .
41 tib a n t r };
42 % We t h e n c h an g e t h e c o m p l i a n c e o f t h e A c h i l l e s t e n do n by c h a n g i n g t h e
11
43 % p a ra m et er Ate ndon o f t h e g a s t r o c n e m i u s and t h e s o l e u s . The d e f a u l t
44 % v a l u e i s 3 5 and a l o w e r v a l u e w i l l r e s u l t i n a more c o m p li a n t t e n d o n .
45 M is c. Atend on = [ 3 5 , 3 5 , 3 5 , 3 5 , 3 5 , 3 5 , 1 5 , 1 5 , 3 5 ] ;
46
47 %% S o l v e t h e p r ob le m
48 switch framework
49 c a s e GPOPS
50 switch formulation actdyn
51 c a s e DeGroote2016
52 switch formulation contdyn
53 c a s e lMtildeState
54 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
55 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
56 SolveMuscleRedundancy lMtildeState GPOPS ( ...
57 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
58 c a s e FtildeState
59 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
60 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
61 SolveMuscleRedundancy FtildeState GPOPS ( . . .
62 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
63 end
64
65 c a s e DeGroote2009
66 switch formulation contdyn
67 c a s e lMtildeState
68 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
69 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
70 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
71 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
72 SolveMuscleRedundancy lMtildeState actdyn GPOPS ( ...
73 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
74 c a s e FtildeState
75 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
76 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
77 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
78 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
79 SolveMuscleRedundancy FtildeState actdyn GPOPS ( ...
80 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
81 end
82 end
83 c a s e CasADi
84 switch formulation actdyn
85 c a s e DeGroote2016
86 switch formulation contdyn
87 c a s e lMtildeState
88 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
89 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
90 S o l v e M u s c l e R e d u n d a n c y l M t i l d e S t a t e C a s A D i ( ...
91 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
92 c a s e FtildeState
93 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
94 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
95 SolveMuscleRedundancy FtildeState CasADi ( ...
96 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
97 end
98
99 c a s e DeGroote2009
100 switch formulation contdyn
101 case lMtildeState
102 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
103 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
104 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
105 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
106 SolveMuscleRedundancy lMtildeState actdyn CasADi ( ...
107 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
108 case FtildeState
109 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
110 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
111 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
112 O p t In f o a c t d yn , D a t S t o r e a c t d y n ]= . . .
113 SolveMuscleRedundancy FtildeState actdyn CasADi ( . . .
114 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
115 end
116 end
117 end
9.4 OpenSim installation example Gait23dof54m
1 c l e a r a l l ; c l o s e a l l ; c l c
2
3%% Choose o p t i m i z a t i o n f ra mewor k
4% fra m ework = GPOPS ;
5 framework = CasADi ;
6
7%% Choose c o n t r a c t i o n dynami cs f o r mu l a t io n
8% f o r m u l a t i o n c o n t d y n = l M t i l d e S t a t e ;
9 f o r m u l a t i o n c o n t d y n = FtildeState ;
10
12
11 %% Choose a c t i v a t i o n dynamics f o r m u l a t io n
12 % f o r m u l a t i o n a c t d y n = D eGr oote 201 6 ;
13 formulation actdyn = DeGroote2009 ;
14
15 %% Example
16 % add main f o l d e r and s u b f o l d e r to matl a b pat h ( i n s t a l l a t i o n )
17 f i l e p a t h=which ( Example Gait23dof54m.m ) ;
18 [ DirExample ,¬,¬]= f i l e p a r t s ( f i l e p a t h ) ;
19 [ Di rEx amp le2 , ¬,¬]= f i l e p a r t s ( DirExample ) ;
20 [ MainD ir , ¬]= f i l e p a r t s ( DirExampl e2 ) ;
21 add path ( ge n p a t h ( MainDir ) ) ;
22
23 % Needed I np u t Arguments
24 Datapath= C : \OpenSim 3 . 3 \Models\G ai t2354 Sim bo dy \OutputReference ;
25 I K p a th= f u l l f i l e ( Datapath , subject01 walk1 ik.mot ) ;
26 I D p ath= f u l l f i l e ( Datapath , R e s u l t s I n v e r s e D y n a m i c s ,inverse dynamics.sto ) ;
27 m od e l pa t h= f u l l f i l e ( Datapath , subject01 scaledOnly.osim ) ;
28 t im e =[0 . 7 1 . 4 ] ; % P art o f t h e r i g h t s t a n c e p ha se
29 Out p ath= f u l l f i l e ( MainDir , Examples ,OpenSimInstallation Gait23dof54m ,Results ) ;
30
31 Misc.DofNames Input={ankle angle r ,k n e e a n g l e r ,h i p f l e x i o n r };
32 Misc.MuscleNames Input={};% S e l e c t s a l l m usc le s f o r i np u t DOFs when empty
33
34 %% S o l v e t h e p r ob le m
35 switch framework
36 c a s e GPOPS
37 switch formulation actdyn
38 c a s e DeGroote2016
39 switch formulation contdyn
40 c a s e lMtildeState
41 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
42 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
43 SolveMuscleRedundancy lMtildeState GPOPS ( ...
44 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
45 c a s e FtildeState
46 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
47 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
48 SolveMuscleRedundancy FtildeState GPOPS ( . . .
49 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
50 end
51
52 c a s e DeGroote2009
53 switch formulation contdyn
54 c a s e lMtildeState
55 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
56 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
57 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
58 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
59 SolveMuscleRedundancy lMtildeState actdyn GPOPS ( ...
60 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
61 c a s e FtildeState
62 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
63 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
64 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
65 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
66 SolveMuscleRedundancy FtildeState actdyn GPOPS ( ...
67 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
68 end
69 end
70 c a s e CasADi
71 switch formulation actdyn
72 c a s e DeGroote2016
73 switch formulation contdyn
74 c a s e lMtildeState
75 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
76 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
77 S o l v e M u s c l e R e d u n d a n c y l M t i l d e S t a t e C a s A D i ( ...
78 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
79 c a s e FtildeState
80 [ Time , M Ex ci ta t io n , M Ac tivat i on , RA c t iv a t io n , T F o r c e t i l d e , ...
81 TFo rc e , l M t i l d e , lM , Mus cleNames , Op tI nfo , D a tS t o re ]= ...
82 SolveMuscleRedundancy FtildeState CasADi ( ...
83 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
84 end
85
86 c a s e DeGroote2009
87 switch formulation contdyn
88 c a s e lMtildeState
89 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
90 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
91 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
92 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
93 SolveMuscleRedundancy lMtildeState actdyn CasADi ( ...
94 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
95 c a s e FtildeState
96 [ Tim e a ctdyn , M E x c i t a ti on a c t d y n , MA c t i va t io n a c t d y n , ...
97 R A ct i va t io n a c t dy n , T F o r c e t i l d e a c t d y n , TF o rc e a ct d yn , ...
98 l M t i l d e a c t d y n , l M ac td yn , M us cl eNam es a ct dy n , . . .
99 O p t I nf o a c t d yn , D a t S t o r e a c t d y n ]= . . .
100 SolveMuscleRedundancy FtildeState actdyn CasADi ( . . .
101 m o d e l p a t h , IK p a th , ID p a th , time , Out path , Misc ) ;
102 end
103 end
104 end
13
