Lecture 9 - Motor Control

Introduction to modeling and simulation of human movement

https://github.com/BMClab/bmc/blob/master/courses/ModSim2018.md

  • In class:
In [115]:
import numpy as np
#import pandas as pd
#import pylab as pl
import matplotlib.pyplot as plt
import math

%matplotlib notebook

Muscle properties

In [116]:
Lslack = .223
Lce_o = .093 #optmal l

Fmax = 3000

Initial conditions

In [117]:
LceNorm = .087/Lce_o
Lm = .31
act = 0
t0 = 0
tf = 2.99
h = 1e-3
In [118]:
t = np.arange(t0,tf,h)
F = np.empty(t.shape)
Fkpe = np.empty(t.shape)
fiberLength = np.empty(t.shape)
tendonLength = np.empty(t.shape)
activation = np.zeros(t.shape)
neural = np.zeros(t.shape)

fiberLength[0] = LceNorm * Lce_o
tendonLength[0] = Lm - fiberLength[0]
In [119]:
def computeTendonForce(LseeNorm, Lslack, Lce_o): 
    '''
    Compute Tendon Force
    
    Imputs: 
    
    LseeNorm - Normalized Tendon length
    
    Lslack - slack length of the tendon (not-normalized)
    
    Lce_o - Optimal length of the fiber 
    
    Output:
    
    FTendonNorm - Normalized tendon force 
    '''
    
    Umax = .04
    if LseeNorm<Lslack/Lce_o: 
        FTendonNorm = 0
    else: 
        FTendonNorm = ((LseeNorm-Lslack/Lce_o)/(Umax*Lslack/Lce_o))**2
        
    return FTendonNorm
In [120]:
def computeParallelElementForce(LceNorm): 
    Umax = 1 
    
    if LceNorm<1: 
        FkpeNorm = 0
    else: 
        FkpeNorm = ((LceNorm-1)/Umax)**2
        
    return FkpeNorm 
In [121]:
def computeForceLengthCurve(LceNorm):

    width = .63
    F0 = max([0, (1-((LceNorm-1)/width)**2)])

    return F0
In [122]:
def computeContractilElementDerivative (act, F0, FCE): 
    
    A = 0.25
    Fm_len = 1.4
    v_max = 10
    
    if  FCE <= act*F0: 
        b = act*F0 + FCE/A
        
    else:
        b = (2 + 2/A)*(act*F0*Fm_len + FCE)/(Fm_len - 1)
       
    LceNormdot = (0.25+0.75*act)*v_max*(FCE-act*F0)/b
    
    return LceNormdot    
In [123]:
def computeContractilElementForce (FTendonNorm, FkpeNorm):
    FCE = FTendonNorm - FkpeNorm
    return FCE
In [124]:
def computeTendonLength(Lm, Lce_o, LceNorm):
    LseeNorm = Lm/Lce_o - LceNorm
    return LseeNorm
In [125]:
def  computeActivation (act, n, dt):
    
    tau = 0.12
    
    dactdt = (n-act)/tau
    act = act + dt*dactdt
    
    return act

Simulation - Parallel

In [126]:
for i in range (len(t)):
    #ramp
    if t[i]<=1:
        Lm = 0.31
    elif t[i]>1 and t[i]<2:
        Lm = .31 - .04*(t[i]-1)
        #print(Lm)
    
    if t[i]<=1 and np.random.rand(1)< 0.3:
        n = 0.5 
    elif t[i]>1.5 and t[i]<2.5 and np.random.rand(1)< 0.5:
        n = 0.5
    elif t[i]>=2.5 and np.random.rand(1)< 0.01:
        n = 0.5
    else:
        n = 0.5 
        
    ##################################################
    LseeNorm = computeTendonLength(Lm, Lce_o, LceNorm)
    
    FTendonNorm = computeTendonForce(LseeNorm, Lslack, Lce_o)    
    
    FkpeNorm = computeParallelElementForce(LceNorm)      
    
    F0 = computeForceLengthCurve (LceNorm)
    
    FCE = computeContractilElementForce (FTendonNorm, FkpeNorm)
    
    act = computeActivation(act, n, h)
          
    LceNormdot = computeContractilElementForce (act, F0, FCE)
    
    
    # --- Euler integration step
    LceNorm = LceNorm + h*LceNormdot
    
    F[i] = FTendonNorm * Fmax
    Fkpe[i] = FkpeNorm * Fmax
    fiberLength[i] = LceNorm * Lce_o
    tendonLength[i] = LseeNorm * Lce_o
    neural[i] = n
    
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-126-8b186d5c87f6> in <module>()
     29     act = computeActivation(act, n, h)
     30 
---> 31     LceNormdot = computeContractilElementForce (act, F0, FCE)
     32 
     33 

TypeError: computeContractilElementForce() takes 2 positional arguments but 3 were given

Plot

In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(6,6), sharex=True)

ax.plot(t,fiberLength, label = 'Fiber')
ax.plot(t,tendonLength, label = 'Tendon')
ax.plot(t,fiberLength+tendonLength, label = 'Fiber+Tendon')
plt.grid()
plt.legend(loc = 'best')
plt.xlabel('time (s)')
plt.ylabel('Length (m)')
plt. tight_layout()