# Task for lecture 20

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.


```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
```

### Muscle force estimation

Consider the following main flexors of the elbow joint (see figure below): biceps long head, biceps short head, and brachialis. Suppose that the elbow net joint torque determined using inverse dynamics is 20 Nm (flexor). How much each of these muscles contributed to the net torque?

<div class='center-align'><figure> <figcaption><center><i>Figure. A view in OpenSim of the arm26 model showing three elbow flexors (Biceps long and short heads and Brachialis).</i></center></figcaption> </figure></div>

For the optimization, we will need experimental data for the moment arm, maximum moment, and *pcsa* of each muscle. Let's import these data from the OpenSim arm26 model:


```python
# time elbow_flexion BIClong BICshort BRA
r_ef = np.loadtxt('./../data/r_elbowflexors.mot', skiprows=7)
f_ef = np.loadtxt('./../data/f_elbowflexors.mot', skiprows=7)
```

The maximum isometric force of these muscles are defined in the arm26 model as: Biceps long head: 624.3 N, Biceps short head: 435.56 N, and Brachialis: 987.26 N. Let's compute the mamimum torques that each muscle could produce considering a static situation at the different elbow flexion angles:


```python
m_ef = r_ef*1
m_ef[:, 2:] = r_ef[:, 2:]*f_ef[:, 2:]
```

These data don't have the *pcsa* value of each muscle. We will estimate the *pcsa* considering that the amount of maximum muscle force generated per area is constant and equal to 50N/cm$^2$. Consequently, the *pcsa* (in cm$^2$) for each muscle is:


```python
a_ef = np.array([624.3, 435.56, 987.26])/50  # 50 N/cm2
print(a_ef)
```

    [12.486   8.7112 19.7452]


### Solution of the optimization problem

For the present case, we are dealing with a problem of minimization, multidimensional (function of several variables), nonlinear, constrained, and we can't assume that the cost function is convex. Numerical optimization is hardly a simple task. There are many different algorithms and public and commercial software for performing optimization. For instance, look at [NEOS Server](http://www.neos-server.org/neos/), a free internet-based service for solving numerical optimization problems.  
We will solve the present problem using the [scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize) package which provides several optimization algorithms. We will use the function `minimize`:

```python
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
"""Minimization of scalar function of one or more variables."""
```

Now, let's write Python functions for each cost function:


```python
from scipy.optimize import minimize
```


```python
def cf_f1(x):
    """Cost function: sum of forces."""  
    return x[0] + x[1] + x[2]

def cf_f2(x):
    """Cost function: sum of forces squared."""
    return x[0]**2 + x[1]**2 + x[2]**2

def cf_fpcsa2(x, a):
    """Cost function: sum of squared muscle stresses."""
    return (x[0]/a[0])**2 + (x[1]/a[1])**2 + (x[2]/a[2])**2

def cf_fmmax3(x, m):
    """Cost function: sum of cubic forces normalized by moments."""
    return (x[0]/m[0])**3 + (x[1]/m[1])**3 + (x[2]/m[2])**3
```

Let's also define the Jacobian for each cost function (which is an optional parameter for the optimization):


```python
def cf_f1d(x):
    """Derivative of cost function: sum of forces."""
    dfdx0 = 1
    dfdx1 = 1
    dfdx2 = 1
    return np.array([dfdx0, dfdx1, dfdx2])

def cf_f2d(x):
    """Derivative of cost function: sum of forces squared."""
    dfdx0 = 2*x[0]
    dfdx1 = 2*x[1]
    dfdx2 = 2*x[2]
    return np.array([dfdx0, dfdx1, dfdx2])

def cf_fpcsa2d(x, a):
    """Derivative of cost function: sum of squared muscle stresses."""
    dfdx0 = 2*x[0]/a[0]**2
    dfdx1 = 2*x[1]/a[1]**2
    dfdx2 = 2*x[2]/a[2]**2
    return np.array([dfdx0, dfdx1, dfdx2])

def cf_fmmax3d(x, m):
    """Derivative of cost function: sum of cubic forces normalized by moments."""
    dfdx0 = 3*x[0]**2/m[0]**3
    dfdx1 = 3*x[1]**2/m[1]**3
    dfdx2 = 3*x[2]**2/m[2]**3
    return np.array([dfdx0, dfdx1, dfdx2])
```

Perform the optimization using the previous initial condition looked in class

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. 


```python
M = 20  # desired torque at the elbow
iang = 69  # which will give the closest value to 90 degrees
#iang = 34
r  = r_ef[iang, 2:]
f0 = f_ef[iang, 2:]
a  = a_ef
m  = m_ef[iang, 2:]
x0 = f_ef[iang, 2:]/10  # far from the correct answer for the sum of torques
print('M =', M)
print('x0 =', x0)
print('r * x0 =', np.sum(r*x0))
```

    M = 20
    x0 = [57.51311369 36.29974032 89.6470056 ]
    r * x0 = 6.622004446068322


Inequality constraints (such as boundaries in our problem) can be entered with the parameter `bounds` to the `minimize` function:


```python
bnds = ((0, f0[0]), (0, f0[1]), (0, f0[2]))
```

Equality constraints (such as the sum of torques should equals the desired torque in our problem), as well as inequality constraints, can be entered with the parameter `constraints` to the `minimize` function (and we can also opt to enter the Jacobian of these constraints):


```python
# use this in combination with the parameter bounds:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)})
```


```python
# to enter everything as constraints:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[0]-x[0],
         'jac' : lambda x, r, f0, M: np.array([-1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[1]-x[1],
         'jac' : lambda x, r, f0, M: np.array([0, -1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[2]-x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, -1]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[0],
         'jac' : lambda x, r, f0, M: np.array([1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[1],
         'jac' : lambda x, r, f0, M: np.array([0, 1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, 1]), 'args': (r, f0, M)})
```

Although more verbose, if all the Jacobians of the constraints are also informed, this alternative seems better than informing bounds for the optimization process (less error in the final result and less iterations).  

Given the characteristics of the problem, if we use the function `minimize` we are limited to the SLSQP (Sequential Least SQuares Programming) solver.  

Finally, let's run the optimization for the four different cost functions and find the optimal muscle forces:


```python
f1r = minimize(fun=cf_f1, x0=x0, args=(), jac=cf_f1d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 409.59266009953336
                Iterations: 7
                Function evaluations: 7
                Gradient evaluations: 7



```python
f2r = minimize(fun=cf_f2, x0=x0, args=(), jac=cf_f2d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 75657.38479127169
                Iterations: 5
                Function evaluations: 6
                Gradient evaluations: 5



```python
fpcsa2r = minimize(fun=cf_fpcsa2, x0=x0, args=(a,), jac=cf_fpcsa2d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 529.9639777695737
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
fmmax3r = minimize(fun=cf_fmmax3, x0=x0, args=(m,), jac=cf_fmmax3d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 1075.1388931664733
                Iterations: 12
                Function evaluations: 12
                Gradient evaluations: 12



```python
# Storing the results of the first optimization using previous initial conditions
opt1 = np.around(np.vstack((f1r.x, f2r.x, fpcsa2r.x, fmmax3r.x)), 1)
```


```python
dat = np.vstack((np.around(r*100,1), np.around(a,1), np.around(f0,0), np.around(m,1)))
er = ['-', '-', '-', '-',
      np.sum(r*f1r.x)-M, np.sum(r*f2r.x)-M, np.sum(r*fpcsa2r.x)-M, np.sum(r*fmmax3r.x)-M]
data = np.vstack((np.vstack((dat, opt1)).T, er)).T

rows = ['$\text{Moment arm}\;[cm]$', '$pcsa\;[cm^2]$', '$F_{max}\;[N]$', '$M_{max}\;[Nm]$',
        '$\sum F_i$', '$\sum F_i^2$', '$\sum(F_i/pcsa_i)^2$', '$\sum(F_i/M_{max,i})^3$']
cols = ['Biceps long head', 'Biceps short head', 'Brachialis', 'Error in M']
df = pd.DataFrame(data, index=rows, columns=cols)
print('\Solving the distribution problem for a desired torque of 20Nm and at 90 elbow flexion')
df
```

    \Solving the distribution problem for a desired torque of 20Nm and at 90 elbow flexion





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Biceps long head</th>
      <th>Biceps short head</th>
      <th>Brachialis</th>
      <th>Error in M</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\text{Moment arm}\;[cm]$</th>
      <td>4.9</td>
      <td>4.9</td>
      <td>2.3</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$pcsa\;[cm^2]$</th>
      <td>12.5</td>
      <td>8.7</td>
      <td>19.7</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$F_{max}\;[N]$</th>
      <td>575.0</td>
      <td>363.0</td>
      <td>896.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$M_{max}\;[Nm]$</th>
      <td>28.1</td>
      <td>17.7</td>
      <td>20.4</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$\sum F_i$</th>
      <td>215.4</td>
      <td>194.2</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum F_i^2$</th>
      <td>184.7</td>
      <td>184.7</td>
      <td>86.1</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/pcsa_i)^2$</th>
      <td>201.7</td>
      <td>98.2</td>
      <td>235.2</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/M_{max,i})^3$</th>
      <td>241.1</td>
      <td>120.9</td>
      <td>102.0</td>
      <td>-3.552713678800501e-15</td>
    </tr>
  </tbody>
</table>
</div>



## 3. 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.  
Set different initial conditions such as consider a desired torque at the elbow of 40 Nm and a slightly higher initial sum of torques.


```python
M = 40  # desired torque at the elbow
iang = 69  # which will give the closest value to 90 degrees
r1  = r_ef[iang, 2:]
f01 = f_ef[iang, 2:]
a1  = a_ef
m1  = m_ef[iang, 2:]
x01 = f_ef[iang, 2:]/5  # far from the correct answer for the sum of torques
print('M =', M)
print('x0 =', x0)
print('r * x0 =', np.sum(r*x0))
```

    M = 40
    x0 = [57.51311369 36.29974032 89.6470056 ]
    r * x0 = 6.622004446068322


Inequality constraints (such as boundaries in our problem) can be entered with the parameter `bounds` to the `minimize` function:


```python
bnds = ((0, f0[0]), (0, f0[1]), (0, f0[2]))
```

Equality constraints (such as the sum of torques should equals the desired torque in our problem), as well as inequality constraints, can be entered with the parameter `constraints` to the `minimize` function (and we can also opt to enter the Jacobian of these constraints):


```python
# use this in combination with the parameter bounds:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)})
```


```python
# to enter everything as constraints:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[0]-x[0],
         'jac' : lambda x, r, f0, M: np.array([-1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[1]-x[1],
         'jac' : lambda x, r, f0, M: np.array([0, -1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[2]-x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, -1]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[0],
         'jac' : lambda x, r, f0, M: np.array([1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[1],
         'jac' : lambda x, r, f0, M: np.array([0, 1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, 1]), 'args': (r, f0, M)})
```

Although more verbose, if all the Jacobians of the constraints are also informed, this alternative seems better than informing bounds for the optimization process (less error in the final result and less iterations).  

Given the characteristics of the problem, if we use the function `minimize` we are limited to the SLSQP (Sequential Least SQuares Programming) solver.  

Finally, let's run the optimization for the four different cost functions and find the optimal muscle forces:


```python
f1r = minimize(fun=cf_f1, x0=x0, args=(), jac=cf_f1d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 819.1853201990796
                Iterations: 10
                Function evaluations: 10
                Gradient evaluations: 10



```python
f2r = minimize(fun=cf_f2, x0=x0, args=(), jac=cf_f2d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 302704.8453888565
                Iterations: 5
                Function evaluations: 6
                Gradient evaluations: 5



```python
fpcsa2r = minimize(fun=cf_fpcsa2, x0=x0, args=(a,), jac=cf_fpcsa2d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 2119.855911078297
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
fmmax3r = minimize(fun=cf_fmmax3, x0=x0, args=(m,), jac=cf_fmmax3d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 8601.111145156132
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
opt2 = np.around(np.vstack((f1r.x, f2r.x, fpcsa2r.x, fmmax3r.x)), 1)
```

Displaying the results for the **second** optimization using different initial conditions basically changing the desired torque (twice as the previous) and initial sum of torques.


```python
dat = np.vstack((np.around(r*100,1), np.around(a,1), np.around(f0,0), np.around(m,1)))
er = ['-', '-', '-', '-',
      np.sum(r*f1r.x)-M, np.sum(r*f2r.x)-M, np.sum(r*fpcsa2r.x)-M, np.sum(r*fmmax3r.x)-M]
data = np.vstack((np.vstack((dat, opt2)).T, er)).T

rows = ['$\text{Moment arm}\;[cm]$', '$pcsa\;[cm^2]$', '$F_{max}\;[N]$', '$M_{max}\;[Nm]$',
        '$\sum F_i$', '$\sum F_i^2$', '$\sum(F_i/pcsa_i)^2$', '$\sum(F_i/M_{max,i})^3$']
cols = ['Biceps long head', 'Biceps short head', 'Brachialis', 'Error in M']
df = pd.DataFrame(data, index=rows, columns=cols)
print('\nSolving the distribution problem for a desired torque of 40Nm')
df
```

    
    Solving the distribution problem for a desired torque of 40Nm





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Biceps long head</th>
      <th>Biceps short head</th>
      <th>Brachialis</th>
      <th>Error in M</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\text{Moment arm}\;[cm]$</th>
      <td>4.9</td>
      <td>4.9</td>
      <td>2.3</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$pcsa\;[cm^2]$</th>
      <td>12.5</td>
      <td>8.7</td>
      <td>19.7</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$F_{max}\;[N]$</th>
      <td>575.0</td>
      <td>363.0</td>
      <td>896.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$M_{max}\;[Nm]$</th>
      <td>28.1</td>
      <td>17.7</td>
      <td>20.4</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$\sum F_i$</th>
      <td>456.2</td>
      <td>363.0</td>
      <td>0.0</td>
      <td>-7.105427357601002e-15</td>
    </tr>
    <tr>
      <th>$\sum F_i^2$</th>
      <td>374.7</td>
      <td>363.0</td>
      <td>174.7</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/pcsa_i)^2$</th>
      <td>403.4</td>
      <td>196.4</td>
      <td>470.5</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/M_{max,i})^3$</th>
      <td>482.2</td>
      <td>241.8</td>
      <td>204.1</td>
      <td>-7.105427357601002e-15</td>
    </tr>
  </tbody>
</table>
</div>



## 3. Regarding the distribution problem for the elbow muscles presented in this text:    
    b. Test other values for the elbow angle where the results are likely to change.
    
Coinsider and angle of 45 degrees.


```python
M = 20  # desired torque at the elbow
iang = 34  # which will give the closest value to 45 degrees
#iang = 34
r  = r_ef[iang, 2:]
f0 = f_ef[iang, 2:]
a  = a_ef
m  = m_ef[iang, 2:]
x0 = f_ef[iang, 2:]/10  # far from the correct answer for the sum of torques
print('M =', M)
print('x0 =', x0)
print('r * x0 =', np.sum(r*x0))
```

    M = 20
    x0 = [62.38587639 43.36183466 98.2329904 ]
    r * x0 = 4.799197847242849


Inequality constraints (such as boundaries in our problem) can be entered with the parameter `bounds` to the `minimize` function:


```python
bnds = ((0, f0[0]), (0, f0[1]), (0, f0[2]))
```

Equality constraints (such as the sum of torques should equals the desired torque in our problem), as well as inequality constraints, can be entered with the parameter `constraints` to the `minimize` function (and we can also opt to enter the Jacobian of these constraints):


```python
# use this in combination with the parameter bounds:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)})
```


```python
# to enter everything as constraints:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[0]-x[0],
         'jac' : lambda x, r, f0, M: np.array([-1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[1]-x[1],
         'jac' : lambda x, r, f0, M: np.array([0, -1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[2]-x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, -1]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[0],
         'jac' : lambda x, r, f0, M: np.array([1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[1],
         'jac' : lambda x, r, f0, M: np.array([0, 1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, 1]), 'args': (r, f0, M)})
```

Although more verbose, if all the Jacobians of the constraints are also informed, this alternative seems better than informing bounds for the optimization process (less error in the final result and less iterations).  

Given the characteristics of the problem, if we use the function `minimize` we are limited to the SLSQP (Sequential Least SQuares Programming) solver.  

Finally, let's run the optimization for the four different cost functions and find the optimal muscle forces:


```python
f1r = minimize(fun=cf_f1, x0=x0, args=(), jac=cf_f1d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 590.5532154884442
                Iterations: 7
                Function evaluations: 7
                Gradient evaluations: 7



```python
f2r = minimize(fun=cf_f2, x0=x0, args=(), jac=cf_f2d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 163425.69129493315
                Iterations: 5
                Function evaluations: 6
                Gradient evaluations: 5



```python
fpcsa2r = minimize(fun=cf_fpcsa2, x0=x0, args=(a,), jac=cf_fpcsa2d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 1227.8544785495337
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
fmmax3r = minimize(fun=cf_fmmax3, x0=x0, args=(m,), jac=cf_fmmax3d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 7770.433925704175
                Iterations: 10
                Function evaluations: 10
                Gradient evaluations: 10



```python
opt3 = np.around(np.vstack((f1r.x, f2r.x, fpcsa2r.x, fmmax3r.x)), 1)
```

**Results of the optimization problem considering a desired torque of 20Nm and 45 degrees elbow flexion**


```python
dat = np.vstack((np.around(r*100,1), np.around(a,1), np.around(f0,0), np.around(m,1)))
er = ['-', '-', '-', '-',
      np.sum(r*f1r.x)-M, np.sum(r*f2r.x)-M, np.sum(r*fpcsa2r.x)-M, np.sum(r*fmmax3r.x)-M]
data = np.vstack((np.vstack((dat, opt3)).T, er)).T

rows = ['$\text{Moment arm}\;[cm]$', '$pcsa\;[cm^2]$', '$F_{max}\;[N]$', '$M_{max}\;[Nm]$',
        '$\sum F_i$', '$\sum F_i^2$', '$\sum(F_i/pcsa_i)^2$', '$\sum(F_i/M_{max,i})^3$']
cols = ['Biceps long head', 'Biceps short head', 'Brachialis', 'Error in M']
df = pd.DataFrame(data, index=rows, columns=cols)
print('\nSolving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion')
df
```

    
    Solving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Biceps long head</th>
      <th>Biceps short head</th>
      <th>Brachialis</th>
      <th>Error in M</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\text{Moment arm}\;[cm]$</th>
      <td>3.4</td>
      <td>3.4</td>
      <td>1.2</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$pcsa\;[cm^2]$</th>
      <td>12.5</td>
      <td>8.7</td>
      <td>19.7</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$F_{max}\;[N]$</th>
      <td>624.0</td>
      <td>434.0</td>
      <td>982.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$M_{max}\;[Nm]$</th>
      <td>21.1</td>
      <td>14.7</td>
      <td>12.2</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$\sum F_i$</th>
      <td>304.8</td>
      <td>285.8</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum F_i^2$</th>
      <td>276.7</td>
      <td>276.7</td>
      <td>101.3</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/pcsa_i)^2$</th>
      <td>324.1</td>
      <td>157.8</td>
      <td>296.7</td>
      <td>-3.552713678800501e-15</td>
    </tr>
    <tr>
      <th>$\sum(F_i/M_{max,i})^3$</th>
      <td>352.3</td>
      <td>204.1</td>
      <td>93.3</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>



## 4. 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.  


```python
M = 10  # desired torque at the elbow
iang = 34  # which will give the closest value to 45 degrees
r  = np.array([0.02, 0.05, 0.04])
f0 = np.array([1000, 250, 700])
a  = np.array([33, 8, 23])
m  = r*f0
x0 = f0/20  # far from the correct answer for the sum of torques
print('M =', M)
print('x0 =', x0)
print('r * x0 =', np.sum(r*x0))
```

    M = 10
    x0 = [50.  12.5 35. ]
    r * x0 = 3.0250000000000004


Inequality constraints (such as boundaries in our problem) can be entered with the parameter `bounds` to the `minimize` function:


```python
bnds = ((0, f0[0]), (0, f0[1]), (0, f0[2]))
```

Equality constraints (such as the sum of torques should equals the desired torque in our problem), as well as inequality constraints, can be entered with the parameter `constraints` to the `minimize` function (and we can also opt to enter the Jacobian of these constraints):


```python
# use this in combination with the parameter bounds:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)})
```


```python
# to enter everything as constraints:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[0]-x[0],
         'jac' : lambda x, r, f0, M: np.array([-1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[1]-x[1],
         'jac' : lambda x, r, f0, M: np.array([0, -1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[2]-x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, -1]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[0],
         'jac' : lambda x, r, f0, M: np.array([1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[1],
         'jac' : lambda x, r, f0, M: np.array([0, 1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, 1]), 'args': (r, f0, M)})
```

Although more verbose, if all the Jacobians of the constraints are also informed, this alternative seems better than informing bounds for the optimization process (less error in the final result and less iterations).  

Given the characteristics of the problem, if we use the function `minimize` we are limited to the SLSQP (Sequential Least SQuares Programming) solver.  

Finally, let's run the optimization for the four different cost functions and find the optimal muscle forces:


```python
f1r = minimize(fun=cf_f1, x0=x0, args=(), jac=cf_f1d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 200.00000000000276
                Iterations: 12
                Function evaluations: 12
                Gradient evaluations: 12



```python
f2r = minimize(fun=cf_f2, x0=x0, args=(), jac=cf_f2d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 22222.22222222222
                Iterations: 5
                Function evaluations: 6
                Gradient evaluations: 5



```python
fpcsa2r = minimize(fun=cf_fpcsa2, x0=x0, args=(a,), jac=cf_fpcsa2d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 69.34812760055478
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
fmmax3r = minimize(fun=cf_fmmax3, x0=x0, args=(m,), jac=cf_fmmax3d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 267.80139010440996
                Iterations: 13
                Function evaluations: 13
                Gradient evaluations: 13



```python
opt4 = np.around(np.vstack((f1r.x, f2r.x, fpcsa2r.x, fmmax3r.x)), 1)
```

**Results of the optimization problem considering a desired torque of 20Nm and 45 degrees elbow flexion**


```python
dat = np.vstack((np.around(r*100,1), np.around(a,1), np.around(f0,0), np.around(m,1)))
er = ['-', '-', '-', '-',
      np.sum(r*f1r.x)-M, np.sum(r*f2r.x)-M, np.sum(r*fpcsa2r.x)-M, np.sum(r*fmmax3r.x)-M]
data = np.vstack((np.vstack((dat, opt4)).T, er)).T

rows = ['$\text{Moment arm}\;[cm]$', '$pcsa\;[cm^2]$', '$F_{max}\;[N]$', '$M_{max}\;[Nm]$',
        '$\sum F_i$', '$\sum F_i^2$', '$\sum(F_i/pcsa_i)^2$', '$\sum(F_i/M_{max,i})^3$']
cols = ['Biceps long head', 'Biceps short head', 'Brachialis', 'Error in M']
df = pd.DataFrame(data, index=rows, columns=cols)
print('\nSolving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion')
df
```

    
    Solving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Biceps long head</th>
      <th>Biceps short head</th>
      <th>Brachialis</th>
      <th>Error in M</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\text{Moment arm}\;[cm]$</th>
      <td>2.0</td>
      <td>5.0</td>
      <td>4.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$pcsa\;[cm^2]$</th>
      <td>33.0</td>
      <td>8.0</td>
      <td>23.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$F_{max}\;[N]$</th>
      <td>1000.0</td>
      <td>250.0</td>
      <td>700.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$M_{max}\;[Nm]$</th>
      <td>20.0</td>
      <td>12.5</td>
      <td>28.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$\sum F_i$</th>
      <td>0.0</td>
      <td>200.0</td>
      <td>0.0</td>
      <td>-1.7763568394002505e-15</td>
    </tr>
    <tr>
      <th>$\sum F_i^2$</th>
      <td>44.4</td>
      <td>111.1</td>
      <td>88.9</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/pcsa_i)^2$</th>
      <td>151.0</td>
      <td>22.2</td>
      <td>146.7</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/M_{max,i})^3$</th>
      <td>65.5</td>
      <td>51.1</td>
      <td>153.3</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>



  ### b. Test the robustness of the results using different initial values for the muscle forces.  
  ### c. Compare the results for different cost functions.


```python
M = 20  # desired torque at the elbow
iang = 34  # which will give the closest value to 45 degrees
r  = np.array([0.02, 0.05, 0.04])
f0 = np.array([1000, 250, 700])
a  = np.array([33, 8, 23])
m  = r*f0
x0 = f0/10  # far from the correct answer for the sum of torques
print('M =', M)
print('x0 =', x0)
print('r * x0 =', np.sum(r*x0))
```

    M = 20
    x0 = [100.  25.  70.]
    r * x0 = 6.050000000000001


Inequality constraints (such as boundaries in our problem) can be entered with the parameter `bounds` to the `minimize` function:


```python
bnds = ((0, f0[0]), (0, f0[1]), (0, f0[2]))
```

Equality constraints (such as the sum of torques should equals the desired torque in our problem), as well as inequality constraints, can be entered with the parameter `constraints` to the `minimize` function (and we can also opt to enter the Jacobian of these constraints):


```python
# use this in combination with the parameter bounds:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)})
```


```python
# to enter everything as constraints:
cons = ({'type': 'eq',
         'fun' : lambda x, r, f0, M: np.array([r[0]*x[0] + r[1]*x[1] + r[2]*x[2] - M]), 
         'jac' : lambda x, r, f0, M: np.array([r[0], r[1], r[2]]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[0]-x[0],
         'jac' : lambda x, r, f0, M: np.array([-1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[1]-x[1],
         'jac' : lambda x, r, f0, M: np.array([0, -1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: f0[2]-x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, -1]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[0],
         'jac' : lambda x, r, f0, M: np.array([1, 0, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[1],
         'jac' : lambda x, r, f0, M: np.array([0, 1, 0]), 'args': (r, f0, M)},
        {'type': 'ineq', 'fun' : lambda x, r, f0, M: x[2],
         'jac' : lambda x, r, f0, M: np.array([0, 0, 1]), 'args': (r, f0, M)})
```

Although more verbose, if all the Jacobians of the constraints are also informed, this alternative seems better than informing bounds for the optimization process (less error in the final result and less iterations).  

Given the characteristics of the problem, if we use the function `minimize` we are limited to the SLSQP (Sequential Least SQuares Programming) solver.  

Finally, let's run the optimization for the four different cost functions and find the optimal muscle forces:


```python
f1r = minimize(fun=cf_f1, x0=x0, args=(), jac=cf_f1d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 437.4999999999541
                Iterations: 8
                Function evaluations: 8
                Gradient evaluations: 8



```python
f2r = minimize(fun=cf_f2, x0=x0, args=(), jac=cf_f2d,
               constraints=cons, method='SLSQP',
               options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 88888.88888888888
                Iterations: 8
                Function evaluations: 10
                Gradient evaluations: 8



```python
fpcsa2r = minimize(fun=cf_fpcsa2, x0=x0, args=(a,), jac=cf_fpcsa2d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 277.3925104022191
                Iterations: 11
                Function evaluations: 11
                Gradient evaluations: 11



```python
fmmax3r = minimize(fun=cf_fmmax3, x0=x0, args=(m,), jac=cf_fmmax3d,
                   constraints=cons, method='SLSQP',
                   options={'disp': True})
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 2142.4111183696405
                Iterations: 12
                Function evaluations: 12
                Gradient evaluations: 12



```python
opt5 = np.around(np.vstack((f1r.x, f2r.x, fpcsa2r.x, fmmax3r.x)), 1)
```

**Results of the optimization problem considering a desired torque of 20Nm and 45 degrees elbow flexion**


```python
dat = np.vstack((np.around(r*100,1), np.around(a,1), np.around(f0,0), np.around(m,1)))
er = ['-', '-', '-', '-',
      np.sum(r*f1r.x)-M, np.sum(r*f2r.x)-M, np.sum(r*fpcsa2r.x)-M, np.sum(r*fmmax3r.x)-M]
data = np.vstack((np.vstack((dat, opt5)).T, er)).T

rows = ['$\text{Moment arm}\;[cm]$', '$pcsa\;[cm^2]$', '$F_{max}\;[N]$', '$M_{max}\;[Nm]$',
        '$\sum F_i$', '$\sum F_i^2$', '$\sum(F_i/pcsa_i)^2$', '$\sum(F_i/M_{max,i})^3$']
cols = ['Biceps long head', 'Biceps short head', 'Brachialis', 'Error in M']
df = pd.DataFrame(data, index=rows, columns=cols)
print('\nSolving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion')
df
```

    
    Solving the distribution problem for a desired torque of 20Nm and at 45 elbow flexion





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Biceps long head</th>
      <th>Biceps short head</th>
      <th>Brachialis</th>
      <th>Error in M</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\text{Moment arm}\;[cm]$</th>
      <td>2.0</td>
      <td>5.0</td>
      <td>4.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$pcsa\;[cm^2]$</th>
      <td>33.0</td>
      <td>8.0</td>
      <td>23.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$F_{max}\;[N]$</th>
      <td>1000.0</td>
      <td>250.0</td>
      <td>700.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$M_{max}\;[Nm]$</th>
      <td>20.0</td>
      <td>12.5</td>
      <td>28.0</td>
      <td>-</td>
    </tr>
    <tr>
      <th>$\sum F_i$</th>
      <td>-0.0</td>
      <td>250.0</td>
      <td>187.5</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum F_i^2$</th>
      <td>88.9</td>
      <td>222.2</td>
      <td>177.8</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/pcsa_i)^2$</th>
      <td>302.1</td>
      <td>44.4</td>
      <td>293.5</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>$\sum(F_i/M_{max,i})^3$</th>
      <td>130.9</td>
      <td>102.3</td>
      <td>306.7</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>



## Comments

The results show that the estimations for the muscle forces depend on the cost function used in the optimization. Which one is correct? This is a difficult question and it's dependent on the goal of the actual task being modeled. Glitsch and Baumann (1997) investigated the effect of different cost functions on the optimization of walking and running and the predicted muscles forces were compared with the electromyographic activity of the corresponding muscles of the lower limb. They found that, among the analyzed cost functions, the minimization of the sum of squared muscle stresses resulted in the best similarity with the actual electromyographic activity.

In general, one should always test different algorithms and different initial values before settling for the solution found. Downey (2011), Kitchin (2013), and Kiusalaas (2013) present more examples on numerical optimization. The [NEOS Guide](http://neos-guide.org/) is a valuable source of information on this topic and [OpenOpt](http://openopt.org/) is a good alternative software for numerical optimization in Python.

## 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.

## References

- Bergstra B, Bengio Y (2012) [Random Search for Hyper-Parameter Optimization](http://jmlr.csail.mit.edu/papers/volume13/bergstra12a/bergstra12a.pdf). Journal of Machine Learning Research, 13, 281-305.  
- Crowninshield RD, Brand RA (1981) [A physiologically based criterion of muscle force prediction in locomotion](http://www.ncbi.nlm.nih.gov/pubmed/7334039). Journal of Biomechanics, 14, 793–801.  
- Downey AB (2014) [Physical Modeling in MATLAB](http://greenteapress.com/wp/physical-modeling-in-matlab-2e/). 2nd edition. Green Tea Press.  
- Herzog W (1987) [Individual muscle force estimations using a non-linear optimal design](http://www.ncbi.nlm.nih.gov/pubmed/3682873). J Neurosci Methods, 21, 167-179.  
- Glitsch U, Baumann W (1997) [The three-dimensional determination of internal loads in the lower extremity](http://www.ncbi.nlm.nih.gov/pubmed/9456380). Journal of Biomechanics, 30, 1123–1131.  
- Kitchin J (2013) [pycse - Python Computations in Science and Engineering](http://kitchingroup.cheme.cmu.edu/pycse/).  
- Kiusalaas (2013) [Numerical methods in engineering with Python 3](http://books.google.com.br/books?id=aJkXoxxoCoUC). 3rd edition. Cambridge University Press.  
- Nigg BM and Herzog W (2006) [Biomechanics of the Musculo-skeletal System](https://books.google.com.br/books?id=hOIeAQAAIAAJ&dq=editions:ISBN0470017678). 3rd Edition. Wiley.  
- Robertson G, Caldwell G, Hamill J, Kamen G (2013) [Research Methods in Biomechanics](http://books.google.com.br/books?id=gRn8AAAAQBAJ). 2nd Edition. Human Kinetics.  
- Seireg A, Arvikar RJ (1973) [A mathematical model for evaluation of forces in lower extremeties of the musculo-skeletal system](http://www.ncbi.nlm.nih.gov/pubmed/4706941). Journal of Biomechanics, 6,  313–322, IN19–IN20, 323–326.

d7cc1695fda5a64e3c4e4a81b255cfb56ab7950c
ipynb
Jupyter Notebook
notebooks/Task for lecture 20.ipynb
regifukuchi/bmc-1
f4418212664758511bb3f4d4ca2318ac48a55e88
MIT
notebooks/Task for lecture 20.ipynb
regifukuchi/bmc-1
f4418212664758511bb3f4d4ca2318ac48a55e88
MIT
notebooks/Task for lecture 20.ipynb
regifukuchi/bmc-1
f4418212664758511bb3f4d4ca2318ac48a55e88
MIT
Qwen/Qwen-72B
1. YES
2. YES