One-legged Cycling Time Trial Example

[moorepants]

Setup

I'm using this conda env opty-sympy-dev-env.yml:

name: opty-sympy-dev
channels:
  - conda-forge
dependencies:
  - coverage
  - cyipopt >=1.1.0
  - cython >=0.29.19
  - matplotlib >=3.2.0
  - mpmath
  - numpy >=1.19.0
  - numpydoc
  - pytest
  - pytest-cov
  - scipy >=1.5.0
  - sphinx
  - sphinx-gallery

after activating I install dev versions of sympy and opty, for example this can work:

python -m pip install --no-deps git+https://github.com/sympy/sympy.git
python -m pip install --no-deps git+https://github.com/csu-hmc/opty.git

or clone them and do an editable install.

TODO

• Choose and set the muscle parameters (can base on numbers from https://github.com/brocksam/muscle-driven-bicycle-paper/blob/main/src/generate_eom.py#L756)
• Read https://docs.sympy.org/dev/tutorials/physics/biomechanics/biomechanics.html carefully and make sure I'm setting muscle properties to reasonable values.
• Set bounds on states and inputs
• Give some kind of rough guesses for some states and inputs
• Try solving the minimum time problem
• Try solving the minimum time + optimal crank & seat height problem
• Reduce the solution time to < 1min (if possible)
• Add text and plots to make the example complete
• Plot muscle forces (may reveal issues)
• Solve with only joint torques
• Check to see if opty actually generates the instance constraints correctly

Current system sketch:

[moorepants]

First solution that shows I'm getting close:

[moorepants]

Closer:

one-legged-time-trial-third-solution.mp4

[moorepants]

Something is fishy with u1: 1) it can't seem to enforce u1(0) = 0 (stationary start) and 2) there is this jump that always appears, blowing the constraint violations:

[moorepants]

Leaving it at this stage:

one-legged-time-trial-sixth-solution.mp4

There is still something weird that causes the crank speed to jump. Maybe it is a weird config of the muscle. Not sure. I could try solving with only joint torques to see if it exposes something else. But overall good progress. I'll come back to it later in the week.

[moorepants]

There is some jiggle, but that may because I cut the solve for the long distance before optimal was found.

[Neville-N]

It did find a solution the first time but the second time running it plateaued and ran into max iteration count

[moorepants]

One thing to do would be to give a much better initial guess. For example we could give alternating square waves for the four muscles that have the likely period of the crank cycle. Right now you can see that I have a np.random.random call that makes the initial guess random and it may not always converge since that changes.

[moorepants]

Feel free to try to fix that by making a square wave guess for the muscle excitations.

[moorepants]

You can also find the max_iter line in the code and bump that up.

[moorepants]

We could make an oval chain ring with three parameters: angle relative to crank, ellipse width/height to find the optimal oval chain ring or something similar. But we can first see if a crank length and/or seat post height as a free parameter will solve.

[moorepants]

@Neville-N can you see if you can run this on google colab?

[Neville-N]

I will start looking at it this afternoon

[moorepants]

I tried solving with the seat height also added as a free optmization variable and it found a solution (87cm for the leg setup in fb9cf4c). Results:

first-optimal-seat-height-solution-87cm.mp4

[moorepants]

diff for making seat height free (not committing):

diff --git a/examples-gallery/plot_one_legged_time_trial.py b/examples-gallery/plot_one_legged_time_trial.py
index 4e37339..5210be4 100644
--- a/examples-gallery/plot_one_legged_time_trial.py
+++ b/examples-gallery/plot_one_legged_time_trial.py
@@ -534,7 +534,7 @@ q_0 = np.array([q1_0, q2_0, q3_0, q4_0])
 # Crank revolutions are proportional to distance traveled so the race distance
 # is defined by number of crank revolutions.
 crank_revs = 4
-samples_per_rev = 40
+samples_per_rev = 120
 num_nodes = crank_revs*samples_per_rev + 1
 
 h = sm.symbols('h', real=True)
@@ -574,8 +574,11 @@ bounds = {
     knee_bot_mus.e: (0.0, 1.0),
     knee_top_mus.e: (0.0, 1.0),
     h: (0.0, 0.1),
+    ls: (0.6, 1.2),
 }
 
+del par_map[ls]
+
 # %
 # Instantiate the Optimal Control Problem
 # =======================================
@@ -591,7 +594,7 @@ problem = Problem(
     bounds=bounds,
 )
 problem.add_option('nlp_scaling_method', 'gradient-based')
-problem.add_option('max_iter', 1000)
+problem.add_option('max_iter', 4000)
 
 initial_guess = np.random.random(problem.num_free)
 
@@ -608,6 +611,7 @@ initial_guess[1*num_nodes:2*num_nodes] = q2_guess
 initial_guess[4*num_nodes:5*num_nodes] = u1_guess
 initial_guess[5*num_nodes:6*num_nodes] = u2_guess
 initial_guess[-4*num_nodes - 1:] = 0.5  # e
+initial_guess[-2] = 0.8
 initial_guess[-1] = 0.01
 
 problem.plot_trajectories(initial_guess)
@@ -616,9 +620,12 @@ problem.plot_trajectories(initial_guess)
 solution, info = problem.solve(initial_guess)
 xs, us, ps, h_val= parse_free(solution, len(state_vars), 4, num_nodes,
                               variable_duration=True)
-print('Optimal value h:', solution[-1])
+print('Optimal value h:', h_val)
+print('Optimal value ls:', ps)
 print(info['status_msg'])
 
+p_vals[14] = ps[0]
+
 # %%
 problem.plot_objective_value()

[moorepants]

This is mergable as a working example that is complete. I will iterate on details in future PRs.