Model testing

The Matlab function das3test.m contains code for testing of the model. As a user of the model, you may want to test the model to verify that it works as it should, or you may want to see the effect of a change in model parameters in a well-defined test rather than in a real-time simulation.

Contents

Visualization test

The model is placed in a state where all joint angles (q) are zero, the forward kinematic model is executed and the resulting pose is visualized. The result should be Figure 1. You can use the zoom and rotate tools at the top of the Figure window for a closer inspection.

At this point, the user is asked to press the Enter key to successively more all joints. This is shown for the first two DOFs in figure 6. The visualization test verifies that the kinematic model is coded correctly.

To execute this test, type:

das3test('vis')
*****************************************************
*                      DAS3MEX                      *
*  (c) 2010-2012 Case Western Reserve University    *
*****************************************************
Initializing...
Visualization test...
Hit ENTER to move SC_y
Hit ENTER to move SC_z
Hit ENTER to move SC_x
Hit ENTER to move AC_y
Hit ENTER to move AC_z
Hit ENTER to move AC_x
Hit ENTER to move GH_y
Hit ENTER to move GH_z
Hit ENTER to move GH_yy
Hit ENTER to move EL_x
Hit ENTER to move PS_y
Figure 1: Result of visualization test.

Speed test

During this test, the MEX function that contains the model dynamics is executed many times with random inputs, and the time per execution is reported on the screen.

It is best to do this when no other process is running on the computer.

To execute this test, type:

das3test('speed')
*****************************************************
*                      DAS3MEX                      *
*  (c) 2010-2012 Case Western Reserve University    *
*****************************************************
Initializing...
Speed test...
Computation time for each model evaluation:   1.66 ms
Computation time with muscle dynamics and no multibody dynamics:        0.13 ms

The second line reports the computation speed that is due to only the muscle and joint models, without multibody dynamics. This shows that the computation time is mostly due to the multibody dynamics.

Moment arm test

This test verifies that the muscle path models are working correctly. Each joint angle is set at the midpoint of its range of motion, all muscle moment arms are calculated, and at each joint, the muscle with the largest positive or negative moment arm is reported.

To execute this test, type:

das3test('moment arms')
*****************************************************
*                      DAS3MEX                      *
*  (c) 2010-2012 Case Western Reserve University    *
*****************************************************
Initializing...
Moment arm test...
DOF   Largest positive moment arm          Largest negative moment arm
----  -----------------------------        -----------------------------------
SC_y              132.1 mm (serr_ant_3      )         -120.3 mm (trap_scap_7     )
SC_z              116.9 mm (trap_scap_2     )         -200.1 mm (lat_dorsi_4     )
SC_x              164.9 mm (serr_ant_1      )          -87.1 mm (rhomboid_2      )
AC_y               58.8 mm (lat_dorsi_1     )         -118.2 mm (lev_scap_2      )
AC_z              148.8 mm (serr_ant_2      )         -163.8 mm (rhomboid_5      )
AC_x               68.7 mm (serr_ant_6      )          -62.5 mm (pect_min_3      )
GH_y               31.8 mm (pect_maj_t_3    )          -17.7 mm (delt_scap_3     )
GH_z               28.3 mm (delt_scap11     )          -73.7 mm (ter_maj_2       )
GH_yy              20.4 mm (subscap_8       )          -14.9 mm (ter_min_3       )
EL_x               81.6 mm (brachiorad_1    )          -19.7 mm (tric_lat_1      )
PS_y                8.7 mm (pron_quad_3     )          -17.3 mm (supinator_1     )
Figure 2: Moment arm matrix generated by das3test('moment arms').

A figure will show the nonzero elements of the moment arm matrix (Figure 2). The moment arm matrix, with the names of the degrees of freedom and muscles, and moment arms in mm, is also written on the file momentarms.csv which can be examined in Excel or other spreadsheet software.

Equilibrium test

This test finds the passive equilibrium of the model. Equilibrium implies that the state derivatives dx/dt are zero. Passive implies that the muscle excitations u are zero.

The solution is found by an iterative solver. When the equilibrium is found, Matlab reports joint variables, muscle variables and glenohumeral stability for this equilibrium.

In addition, the linearized dynamics are examined at equilibrium, by eigenvalue analysis. The magnitude of the eigenvalues tells us how fast the dynamics are. If there are eigenvalues with a magnitude of $10^6$, which corresponds to a microsecond time scale, this means that the differential equation is stiff and requires implicit methods for simulation. The real parts of all eigenvalues have to be negative for the equilibrium to be stable.

To execute this test, type:

das3test('equilibrium')
*****************************************************
*                      DAS3MEX                      *
*  (c) 2010-2012 Case Western Reserve University    *
*****************************************************
Initializing...
Equilibrium test...
Equilibrium eval # 1: Norm of f:  1.79963e-12
Equilibrium eval # 2: Norm of f:  3.30170e-12
Equilibrium eval # 3: Norm of f:  3.30164e-12
Equilibrium eval # 4: Norm of f:  4.26102e-12
Equilibrium eval # 5: Norm of f:  4.26102e-12
Equilibrium eval # 6: Norm of f:  4.26108e-12
Equilibrium eval # 7: Norm of f:  4.26401e-12
Equilibrium eval # 8: Norm of f:  4.26401e-12
Equilibrium eval # 9: Norm of f:  4.26401e-12
Equilibrium eval # 10: Norm of f:  2.61000e-12
Equilibrium eval # 11: Norm of f:  2.62967e-12
Equilibrium eval # 12: Norm of f:  2.58002e-12
Equilibrium eval # 13: Norm of f:  2.59973e-12
Equilibrium eval # 14: Norm of f:  2.59973e-12
Equilibrium eval # 15: Norm of f:  3.35742e-12
Equilibrium eval # 16: Norm of f:  3.83374e-12
Equilibrium eval # 17: Norm of f:  3.91520e-12
Equilibrium eval # 18: Norm of f:  3.07880e-12
Equilibrium eval # 19: Norm of f:  3.07880e-12
Equilibrium eval # 20: Norm of f:  3.07880e-12
Equilibrium eval # 21: Norm of f:  3.07880e-12
Equilibrium eval # 22: Norm of f:  3.07880e-12
Equilibrium eval # 23: Norm of f:  3.07880e-12
Equilibrium eval # 24: Norm of f:  3.07880e-12
Equilibrium eval # 25: Norm of f:  3.07880e-12
Equilibrium eval # 26: Norm of f:  3.07880e-12
Equilibrium eval # 27: Norm of f:  3.07880e-12
Equilibrium eval # 28: Norm of f:  3.07880e-12
Equilibrium eval # 29: Norm of f:  1.79963e-12
Norm of f was: 1.7996e-12
This must be close to zero, otherwise it is not an equilibrium
Equilibrium state was stored on equilibrium.mat
Equilibrium posture was stored on eq_position.sto
 
-------------Equilibrium state-------------
 
DOF               angle(deg)   ang.vel(deg/s)   moment(Nm)  
--------------- -------------- -------------- --------------
SC_y              -21.246         -0.000          4.043
SC_z                7.551         -0.000        -12.569
SC_x               12.761         -0.000         26.385
AC_y               35.555         -0.000        -12.303
AC_z              -24.940         -0.000         40.659
AC_x              -17.947         -0.000         18.458
GH_y               24.878          0.000         -0.044
GH_z               28.549          0.000          1.374
GH_yy              15.642         -0.000         -0.219
EL_x               37.638         -0.000          0.651
PS_y               81.441         -0.000          0.040
 
---------Muscle variables at equilibrium---------
 
Muscle           Lce/Lceopt  Active state   Force(N)  
--------------- ------------ ------------ ------------
trap_scap_1         1.224       -0.000       19.122
trap_scap_2         1.120        0.000        7.820
trap_scap_3         0.855       -0.000       -0.155
trap_scap_4         0.813       -0.000       -0.178
trap_scap_5         0.863       -0.000       -0.127
trap_scap_6         0.942       -0.000       -0.052
trap_scap_7         1.003        0.000        0.005
trap_scap_8         1.165       -0.000       12.104
trap_scap_9         1.270       -0.000       28.915
trap_scap10         1.253       -0.000       16.842
trap_scap11         1.372       -0.000       41.133
trap_clav_1         0.894        0.000       -0.165
trap_clav_2         0.844       -0.000       -0.262
lev_scap_1          0.860       -0.000       -0.221
lev_scap_2          0.801       -0.000       -0.259
pect_min_1          0.834       -0.000       -0.127
pect_min_2          0.869       -0.000       -0.143
pect_min_3          0.805       -0.000       -0.210
pect_min_4          0.849       -0.000       -0.208
rhomboid_1          0.807       -0.000       -0.191
rhomboid_2          0.878       -0.000       -0.131
rhomboid_3          0.757       -0.000       -0.306
rhomboid_4          0.665       -0.000       -0.424
rhomboid_5          0.590        0.000       -0.480
serr_ant_1          1.106        0.000        2.373
serr_ant_2          1.182       -0.000       11.894
serr_ant_3          1.258        0.000       22.699
serr_ant_4          1.479        0.000       81.191
serr_ant_5          1.415        0.000       35.743
serr_ant_6          1.252        0.000       20.002
serr_ant_7          1.012        0.000        0.035
serr_ant_8          1.110        0.000        3.028
serr_ant_9          0.969       -0.000       -0.025
serr_ant10          0.835       -0.000       -0.109
serr_ant11          0.775       -0.000       -0.199
serr_ant12          0.695       -0.000       -0.289
delt_scap_1         1.017        0.000        0.417
delt_scap_2         0.969       -0.000       -0.026
delt_scap_3         0.986       -0.000       -0.009
delt_scap_4         1.076        0.000        2.382
delt_scap_5         1.071        0.000        2.183
delt_scap_6         1.013       -0.000        0.276
delt_scap_7         1.035        0.000        0.746
delt_scap_8         1.007        0.000        0.082
delt_scap_9         0.994       -0.000       -0.002
delt_scap10         0.930       -0.000       -0.049
delt_scap11         0.961       -0.000       -0.026
delt_clav_1         1.350        0.000       51.119
delt_clav_2         0.730       -0.000       -0.247
delt_clav_3         0.745       -0.000       -0.204
delt_clav_4         1.156        0.000        9.188
coracobr_1          1.013       -0.000        0.111
coracobr_2          1.017        0.000        0.153
coracobr_3          1.015       -0.000        0.098
infra_1             0.932       -0.000       -0.043
infra_2             0.978        0.000       -0.016
infra_3             0.993       -0.000       -0.004
infra_4             1.042        0.000       -0.040
infra_5             1.138       -0.000       -0.043
infra_6             1.025       -0.000       -0.056
ter_min_1           1.225       -0.000       -0.037
ter_min_2           1.327       -0.000       -0.040
ter_min_3           1.326        0.000       -0.037
ter_maj_1           0.904       -0.000       -0.135
ter_maj_2           0.918       -0.000       -0.131
ter_maj_3           0.895        0.000       -0.160
ter_maj_4           0.833       -0.000       -0.182
supra_1             1.228        0.000        6.371
supra_2             1.093        0.000        3.139
supra_3             0.972       -0.000       -0.016
supra_4             0.972        0.000       -0.016
subscap_1           0.819       -0.000       -0.100
subscap_2           0.858       -0.000       -0.092
subscap_3           0.814        0.000       -0.122
subscap_4           0.843        0.000       -0.113
subscap_5           0.773        0.000       -0.181
subscap_6           0.868       -0.000       -0.101
subscap_7           0.797        0.000       -0.166
subscap_8           0.786       -0.000       -0.164
subscap_9           0.825       -0.000       -0.126
subscap10           0.826        0.000       -0.118
subscap11           0.742       -0.000       -0.166
bic_l               1.147       -0.000       -0.197
bic_b_1             1.356       -0.000       -0.168
bic_b_2             1.411       -0.000       -0.217
tric_long_1         1.323       -0.000       -0.166
tric_long_2         1.344       -0.000       -0.126
tric_long_3         1.272       -0.000       -0.136
tric_long_4         1.130       -0.000       -0.086
lat_dorsi_1         0.722        0.000       -0.651
lat_dorsi_2         0.756       -0.000       -0.698
lat_dorsi_3         0.752       -0.000       -0.731
lat_dorsi_4         0.791        0.000       -0.670
lat_dorsi_5         0.822       -0.000       -0.556
lat_dorsi_6         0.749       -0.000       -0.853
pect_maj_t_1        0.849       -0.000       -0.285
pect_maj_t_2        0.914       -0.000       -0.158
pect_maj_t_3        0.895       -0.000       -0.194
pect_maj_t_4        0.951       -0.000       -0.090
pect_maj_t_5        1.137       -0.000        8.197
pect_maj_t_6        0.993       -0.000       -0.010
pect_maj_c_1        1.093       -0.000        5.039
pect_maj_c_2        0.977       -0.000       -0.025
tric_med_1          1.110       -0.000       -0.119
tric_med_2          1.146       -0.000       -0.119
tric_med_3          1.223       -0.000       -0.106
tric_med_4          1.130       -0.000       -0.120
tric_med_5          1.109        0.000       -0.073
brachialis_1        1.152        0.000       -0.042
brachialis_2        1.175        0.000       -0.024
brachialis_3        1.102        0.000        0.003
brachialis_4        1.363       -0.000       -0.028
brachialis_5        1.287        0.000       -0.008
brachialis_6        1.149        0.000       -0.024
brachialis_7        1.279        0.000       -0.013
brachiorad_1        1.251       -0.000        0.332
brachiorad_2        1.199       -0.000       -0.001
brachiorad_3        1.367        0.000        0.826
pron_teres_1        1.281        0.000       -0.012
pron_teres_2        1.062        0.000        1.726
supinator_1         0.847       -0.000       -0.041
supinator_2         0.944       -0.000       -0.066
supinator_3         1.158       -0.000       -0.036
supinator_4         1.103       -0.000       -0.025
supinator_5         1.185       -0.000       -0.031
pron_quad_1         1.192       -0.000        1.991
pron_quad_2         1.037       -0.000        0.807
pron_quad_3         0.993        0.000       -0.003
tric_lat_1          1.073       -0.000       -0.096
tric_lat_2          1.042       -0.000       -0.118
tric_lat_3          1.031       -0.000       -0.100
tric_lat_4          0.999       -0.000       -0.117
tric_lat_5          0.849        0.000       -0.062
anconeus_1          0.863        0.000       -0.035
anconeus_2          0.879       -0.000       -0.016
anconeus_3          1.019       -0.000        0.032
anconeus_4          0.891       -0.000       -0.014
anconeus_5          0.785       -0.000       -0.029
 
-----------Other variables---------------
 
Force in glenohumeral joint:
  -18.3646   29.2334   18.0681

Glenohumeral stability value:
   -0.2234

Contact force in Trigonum Spinae of the scapula:
   83.2593  123.6619  217.9835

Contact force in Inferior Angle of the scapula:
   42.9447   28.3786  211.3670

Thorax surface equation solved for Trigonum Spinae of the scapula:
   -0.0231

Thorax surface equation solved for Inferior Angle of the scapula:
   -0.0132

 
----------------Stability analysis----------------
 
Eigenvalues with largest real part (s^-1):
(must all be negative for stability)
   -1.0838   -1.1769   -1.3227   -1.3227   -1.9228

Smallest time constants (s):
1.8509e-06  1.9748e-06  2.3240e-06  3.2712e-06  3.4194e-06
 
The imaginary parts, representing (damped) oscillatory behavior,
have a maximum of ~7374, which is equivalent to a frequency (f=omega/2pi) of  1173.6 Hz.
Figure 3: Passive equilibrium posture and eigenvalues of the dynamic model. Note the logarithmic real axis.

Note that this is for a passive equilibrium state. The eigenvalues will probably be larger when muscles are generating force and tendons contribute to joint stiffness.

Simulation test with real-time solver

The system is simulated with the first-order Rosenbrock solver. This is a fixed-step solver that is suitable for real-time.

At t=0, the model is in its equilibrium state, which is obtained from the file equilibrium.mat. It is therefore important that you have first done the equilibrium test with the same model. The muscle excitations u are then ramped up to reach a maximum (1.0) at t=0.5.

To execute this test, type:

>> das3test(‘simulate’, stepsize)

If stepsize is omitted, a stepsize of 0.003 seconds is used. When the simulation is completed, it displays the result and writes the result in files simulation.mat (full state x as function of time) and simulation.sto (skeleton pose q as function of time, for animation with Opensim).

The results of this simulation are compared to the results from the ode15s test. It is therefore important that test(‘ode15s’) has already been completed with the same model. (The results of ode15s are included in the release: ode15s_simulation.mat (full state x as function of time) and ode15s_simulation.sto (q as function of time, to visualize simulated movement in OpenSim). Differences are reported as RMS error in each joint angle.

das3test('simulate')
*****************************************************
*                      DAS3MEX                      *
*  (c) 2010-2012 Case Western Reserve University    *
*****************************************************
Initializing...
Simulating arm movements...
CPU time per time step:    2.278 ms
Simulation speed is    1.314 times faster than real time
Simulation result has been saved in simulation.mat
Simulated motion has been saved in simulation.sto
RMS errors in each joint angle (deg):
SC_y  :    0.064
SC_z  :    0.042
SC_x  :    0.013
AC_y  :    0.006
AC_z  :    0.016
AC_x  :    0.011
GH_y  :    0.110
GH_z  :    0.057
GH_yy :    0.038
EL_x  :    0.029
PS_y  :    0.140

Overall RMS error:    0.063 degrees

ans =

    0.0629

The rest of the tests are quite slow (an hour or more) so please allow enough time.

Simulation test with Matlab solver

The system can be simulated with a Matlab differential equation solver. This is slow, but we know it is accurate because the Matlab solver will control its time steps to remain below the specified error tolerance. Because it is a stiff system, we use the ode15s solver which is an implicit variable step method.

At t=0, the model is in its equilibrium state, which is obtained from the file equilibrium.mat. It is therefore important that you have first done the equilibrium test with the same model. The muscle excitations u are then ramped up to reach a maximum (1.0) at t=0.5.

To execute this test, type:

>> das3test(‘ode15s’)

At each time step of the simulation, the state derivative $\dot \textbf{x}$ will be solved numerically, which is time consuming. Due to the variable step, computation time is unpredictable (this is why this solver is not suitable for real time) and can be several minutes up to an hour. At regular intervals, the test will report its progress, and some verification that the condition f=0 is being maintained:

Step 784511: Neval = 1202744 -- t =      2.8126797952592 -- Norm(f) = 2.197e-016

The simulated movements are reported as graphs (Figure 4) and also saved in files ode15s_simulation.mat (full state x as function of time) and ode15s_simulation.sto (q as function of time, to visualize simulated movement in OpenSim).

Figure 4: Result of movement simulation test.

To see the simulated movement in OpenSim: (1) Start OpenSim, (2) File->Open Model..., select the model das3.osim, (3) File->Load Motion..., select the motion file ode ode15s_simulation.sto. Then use the playback controls in top right of window to show the movement.

Isometric muscle tests

This test produces isometric moment-angle curves (figure 5) for a joint, one for each value in a range of angles in a second joint.

To execute this test, type:

>> das3test('isometric')

Figure 5: Isometric moment-angle curves for forearm
pronation/supination, one for each value in a range of angles of elbow
flexion/extension

You can also produce one isometric moment-angle curve (figure 6) for each muscle that crosses a joint.

To execute this test, type:

>> das3test('isometric1')

Figure 6: Isometric moment-angle curve for all muscles crossing ACy

Isokinetic muscle test

This test produces isokinetic moment-angular velocity curves for a joint.

To execute this test, type:

>> das3test('isokinetic')

Figure 7: Isokinetic moment-angular velocity curve for SCz

Previous: Model reference
Home: Main