Computed Muscle Control Algorithm

[Takuma Inai]

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

[Dimitar Stanev]

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.

[Takuma Inai]

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