Closed moorepants closed 2 months ago
First solution that shows I'm getting close:
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:
Leaving it at this stage:
https://github.com/csu-hmc/opty/assets/276007/e4ab11a0-6e24-495b-a931-67cda01a4611
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.
I just setup everything and have the code running. Is there anything I should first look at?
I now get an optimal solution. I plotted the muscle forces and saw I had the max forces set too low (which were just guesses). Now I can solve from a standstill up to a good cadence. Here is an example solution:
https://github.com/csu-hmc/opty/assets/276007/427bdd18-f42b-4c62-a9e1-919635538bb9
It now works! Run it and see if it solves for you.
There is some jiggle, but that may because I cut the solve for the long distance before optimal was found.
It did find a solution the first time but the second time running it plateaued and ran into max iteration count
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.
Feel free to try to fix that by making a square wave guess for the muscle excitations.
You can also find the max_iter
line in the code and bump that up.
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.
@Neville-N can you see if you can run this on google colab?
I will start looking at it this afternoon
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:
https://github.com/csu-hmc/opty/assets/276007/d60b8f8b-03f0-42f8-842e-33adec787bb5
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()
This is mergable as a working example that is complete. I will iterate on details in future PRs.
Setup
I'm using this conda env
opty-sympy-dev-env.yml
:after activating I install dev versions of sympy and opty, for example this can work:
or clone them and do an editable install.
TODO
Solve with only joint torquesCheck to see if opty actually generates the instance constraints correctlyCurrent system sketch: