## Task 19 - Motor Control
### Introduction to modeling and simulation of human movement
https://github.com/BMClab/bmc/blob/master/courses/ModSim2018.md

Desiree Miraldo

* Task (for Lecture 19):

Solve problemas 1 and 2 of the notebook [Optimization (Marcos Duarte)](http://nbviewer.jupyter.org/github/BMClab/bmc/blob/master/notebooks/Optimization.ipynb).


```python
# import Python libraries
import numpy as np
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import sympy as sym
from sympy.plotting import plot
import pandas as pd
from IPython.display import display
from IPython.core.display import Math
```

## Exercises

1. Find the extrema in the function $f(x)=x^3-7.5x^2+18x-10$ analytically and determine if they are minimum or maximum.  
2. Find the minimum in the $f(x)=x^3-7.5x^2+18x-10$ using the gradient descent algorithm.  
2. Regarding the distribution problem for the elbow muscles presented in this text:  
    a. Test different initial values for the optimization.  
    b. Test other values for the elbow angle where the results are likely to change.   
    
3. In an experiment to estimate forces of the elbow flexors, through inverse dynamics it was found an elbow flexor moment of 10 Nm.  
Consider the following data for maximum force (F0), moment arm (r), and pcsa (A) of the brachialis, brachioradialis, and biceps brachii muscles: F0 (N): 1000, 250, 700; r (cm): 2, 5, 4; A (cm$^2$): 33, 8, 23, respectively (data from Robertson et al. (2013)).  
    a. Use static optimization to estimate the muscle forces.  
    b. Test the robustness of the results using different initial values for the muscle forces.  
    c. Compare the results for different cost functions.

#### Exercises 1: Find the extrema in the function $f(x)=x^3-7.5x^2+18x-10$ analytically and determine if they are minimum or maximum.  


```python
x = sym.symbols('x')
f = x**3 - 7.5*x**2 + 18*x - 10
fdiff = sym.expand(sym.diff(f, x))
roots = sym.solve(fdiff, x)
display(Math(sym.latex('Roots:') + sym.latex(roots)))
```


$$Roots:\left [ 2.0, \quad 3.0\right ]$$



```python
f  = lambda x: x**3 - 7.5*x**2 + 18*x - 10

print('f(x = {:.1f}) = {:.2f}'.format(roots[0],f(roots[0])))
print('f(x = {:.1f}) = {:.2f}'.format(roots[1],f(roots[1])))
```

    f(x = 2.0) = 4.00
    f(x = 3.0) = 3.50


#### Exercises 2: Find the minimum in the $f(x)=x^3-7.5x^2+18x-10$ using the gradient descent algorithm.  


```python
# From https://en.wikipedia.org/wiki/Gradient_descent
# The local minimum of $f(x)=x^3-7.5x^2+18x-10$

cur_x = 6               # The algorithm starts at x=6
gamma = 0.01            # step size multiplier
precision = 0.00001
step_size = 1           # initial step size
max_iters = 10000       # maximum number of iterations
iters = 0               # iteration counter

f  = lambda x: x**3 - 7.5*x**2 + 18*x - 10  # lambda function for f(x)
df = lambda x: 3*x**2 - 15*x + 18    # lambda function for the gradient of f(x)

while (step_size > precision) and (iters < max_iters):
    prev_x = cur_x
    cur_x -= gamma*df(prev_x)
    step_size = abs(cur_x - prev_x)
    iters+=1

print('True local minimum at {} with function value {}.'.format(3, f(3)))
print('Local minimum by gradient descent at {} with function value {}.'.format(cur_x, f(cur_x)))

from scipy.optimize import minimize
import numpy as np

minimize(f, 6, args=(), method='CG', jac=None, tol=precision, callback=None,
         options={'gtol': 1e-05, 'norm': np.inf, 'eps': 1.4901161193847656e-08,
                  'maxiter': max_iters, 'disp': True, 'return_all': False})
```

    True local minimum at 3 with function value 3.5.
    Local minimum by gradient descent at 3.000323195755751 with function value 3.5000001567170003.
    Optimization terminated successfully.
             Current function value: -4150663493044532703920128.000000
             Iterations: 3
             Function evaluations: 93
             Gradient evaluations: 31





         fun: array([ -4.15066349e+24])
         jac: array([ 0.])
     message: 'Optimization terminated successfully.'
        nfev: 93
         nit: 3
        njev: 31
      status: 0
     success: True
           x: array([ -1.60708619e+08])




```python

```

1bdba485da21c41e85273c2c6ee3708f21f2da79
ipynb
Jupyter Notebook
courses/modsim2018/tasks/Desiree/Task19_MotorControl.ipynb
desireemiraldo/bmc
5e50d806ea2e2d7494035ba120c3cd3a620a156a
MIT
courses/modsim2018/tasks/Desiree/Task19_MotorControl.ipynb
desireemiraldo/bmc
5e50d806ea2e2d7494035ba120c3cd3a620a156a
MIT
courses/modsim2018/tasks/Desiree/Task19_MotorControl.ipynb
desireemiraldo/bmc
5e50d806ea2e2d7494035ba120c3cd3a620a156a
MIT
Qwen/Qwen-72B
1. YES
2. YES