simtk.org

Hi,

I want to completely understand the computed muscle control algorithm. To understand the algorithm, I have been reading a previous study by Thelen et al. (2003).

Generating dynamic simulation of movement using computed muscle control. Journal of Biomechanics. 2003.

Would you tell me whether the algorithm I am imaging is correct?

Note:
"_d1" means first-order time differentiation.
"_d2" means second-order time differentiation.
"[ ]" means size of matrix.

Assumption:
- Number of degrees of freedom => 9
- Number of muscles => 16
Reference for model:
Joint contact forces can be reduced by improving joint moment symmetry. Gait & Posture. 2016.

Preparation (Initial values):
lm(t) [16 x 1] : muscle fiber length
a(t) [16 x 1]: muscle activation
q_sim_d1(t) [9 x 1]: generalized velocity (simulation)
q_sim(t) [9 x 1]: generalized coordinate (simulation)

Assumption.
Please neglect pennation angles.

Algorithm:
Stage 1
q_desired_d2(t+T) = q_exp_d2(t+T) + kv*(q_exp_d1(t) - q_sim_d1(t)) + kp*(q_exp(t) - (q_sim(t)))
Note: kv = 20 and kp = 100

From above equation, calculate q_desired_d2(t+T) [9 x 1].

A*q_desired_d2(t+T) = G(q_sim(t)) + C(q_sim(t), q_sim_d1(t)) + E(q_sim(t),q_sim_d1(t)) + Q'

Note: Q' is generalized forces [9 x 1].
A: Mass system matrix
G: Gravity matrix
C: Coriolis and centripetal forces matrix
E: Matrix of external forces
Q': Generalized forces (e.g., joint torque)

From above equation, calculate Q'.

Next, from generalized forces (Q'), calculate Residual forces/torque (R) [3 x 1] and Joint torque (Tau) [6 x 1].

Stage 2

lm_d1(t) = fv^-1(lm(t), lmt(t), a(t))

Note: lm(t) is initial values.
From above equation, calculate lm_d1(t) [16 x 1].

For each muscle,
f# = a# .* flv(lm#, lm_d1#) + fpassive(lm#)-------Eq. A

Note:
f#: muscle force [16 x 1]
a#: muscle activation [16 x 1]
flv: force-length-velocity relationship

lm# <= lm(t) [16 x 1] i.e., initial value.
lm_d1# is calculated from fv^-1(lm(t), lmt(t), a(t))
lm(t) is initial value [16 x 1], lmt(t) is calculated using joint angles (i.e., q_sim(t)), and a(t) is initial value [16 x 1].

Using Joint torque (Tau) [6 x 1], Eq. A, and static optimization, estimate a#.

Stage 3

u(t) = a# + ku*(a# - a(t))

Note:
a# [16 x 1].
a(t) is initial value [16 x 1].
ku is 10.

From above equation, calculate u(t) [16 x 1].

Stage 4
a_d1 = (u-a)*{(u/tact) + (1-u)/tdeact}, u >= a----------Eq. B
a_d1 = (u-a)/tdeact, u < a----------Eq. B

Using Runge-Kutta method with above equations (i.e., a_d1(u,a)) with u(t) [16 x 1] and a(t) [16 x 1], calculate a(t+T) [16 x 1].

Note: tact is 15 ms and tdeact is 50 ms.

After that, using below formula,

f(t+T) = a(t+T) * flv(lm(t+T), lm_d1(t+T)) + fpassive(lm(t+T))-------Eq. A

with below equation of motion,

q_sim_d2(t+T) = A^-1*{G(q_sim(t+T)) + C(q_sim(t+T), q_sim_d1(t+T)) + E(q_sim(t+T),q_sim_d1(t+T)) + Q(t+T)}

calculate q_sim_d2(t+T).

After that, calculate q_sim(t+2*T) and q_sim_d1(t+2*T) using Runge-Kutta method.

Note: Q(t+T) includes f(t+T) for all muscles.

Could you give some advice please?

Sincerely

Takuma Inai
Institute for Human Movement and Medical Sciences
Niigata University of Health and Welfare

It is best if you could dig into the source code:

https://github.com/opensim-org/opensim- ... l.cpp#L384

Theory and implementation may different in practice.

Thank you so much.
I see the code, and I learn the argorithm from the code.