Matlab version of example2DWalking

[Aaron Fox]

Hi Ross,

From what I've gathered about how models work with most Moco processes is that the whole model gets copied into the Moco study - i.e. if you print out the Moco study or tracking object as a .omoco file you will see that the model (and any changes) have been dumped into there. With the model processor, I think instead of the model getting placed in the study the processor does, and any of the changes that are scheduled to be 'processed' don't happen until you run the study. So if you set-up changes in a model processor, they won't occur to the model until the study is run. You have 'processed' it in your code but perhaps these changes won't occur with this? I would have expected them to though.

Not too sure about your first question, but I know that default properties don't get printed when you output a model from Matlab to it's .osim file. Is it possible that those length approximations are the default for the muscle?

Aaron

[Ross Miller]

Thanks Aaron.

Recently I removed the length_approximation's from the model and the optimization still works well. I guess this would probably only be a problem if the muscle_length = f(Q) function is discontinuous.

Ross

[Ross Miller]

Hi all,

Bumping this thread to provide a link to the SimTk project I've been working on lately for optimal control simulations of gait with Moco. It's in a state where I think it's working reasonably well enough for sharing with others:

https://simtk.org/projects/umocod

The link has a download with the Matlab script, .osim model, and necessary data for performing tracking simulations of walking, and could be easily extended to performing predictive simulations. The recent updated release of Moco (0.4.0) is required. The model is a 2D version of the Rajagopal et al. (2016) model. The work is essentially a more complicated version of the 2D walking example included with Moco:

- Enabled muscle contractile dynamics for all muscles
- Added an articulating toes segment and muscles that control it
- Added additional ground contact elements
- Replaced the lumbar actuator with trunk flexor/extensor muscles
- Adjusted some muscle model parameters with respect to referenced human joint torque-angle data
- Defined muscle-specific activation time constants using muscle mass and fiber type data
- Added some coordinate limit forces on the joints
- Set the code up to track mean human experimental gait data (these data are included)
- Set the code up to simulate a full stride
- Added a periodicity goal for the muscle excitations

I thought it would serve as a nice starting point for users looking to do simulations with a more complicated model. Personally I'm endeavoring to extend it to 3D for a DoD project on prostheses (hence the simulation of a full stride, potentially asymmetric) and would like to add a Moco goal that includes muscle metabolic energy expenditure (maybe this is already in the works?). This version of the model is still reasonably fast, most things I've tried with it have converged in an "overnight" timeframe or less on my modest personal computer (2019 MacBook Pro 16", 2.6-GHz 6-core Intel Core i7 CPU). A 3D version I think will probably require a workstation, or a charred husk of a laptop.

Some known issues:

- Muscle wrapping doesn't work. When any wrap object is tagged under the muscles in the model's XML file, Matlab crashes before IPOPT starts iterating. Chris thinks this might be a Mac-specific bug. Curiously, it works on my Mac when the IPOPT parallelism option is disabled in Moco. The download includes a version of the model with wrapping objects tagged under the muscles. If anyone running Windows or Linux is able to get IPOPT to iterate with that model, please let me know. Lack of wrapping causes questionable muscle kinematics in some cases, although they mostly occur outside the ranges of motion for walking.

- In the simulation result provided (tracking problem), rectus femoris produces a passive force of ~20-40% max isometric force during swing. The Rajagopal model is quite strong so this ends up being a pretty large force in an absolute sense. Not sure if this is necessarily unrealistic but thought it was worth pointing out.

- The implementation of passive joint moments is pretty hand-wavy right now. <passive_fiber_strain_at_one_norm_force> is set to 1.0 for all muscles then <CoordinateLimitForce> restricts the joints to arbitrary end ranges of motion. I'd like to make this less arbitrary, e.g. implement something like Riener & Edrich (1999) [see my recent post on passive joint moments on the main OpenSim forum].

Would appreciate any suggestions etc. on these issues or other issues anyone runs across!

Ross

[Mritula Chandrasekaran]

Hi Brain,
Thanks for the wonderful code.
Can you please clarify regarding the coordinates.sto data. How did you choose or derive it or get it? For predictive simulation, my understanding is that we need the model and specific parameters for simulation like time period, gait speed, etc.

Thanks,
Mritula

[Christopher Dembia]

Nick just told me of all the discussion on this thread. Thank you, everyone, for helping each other Thank you, Brian, for mentioning that old forum thread on implicit equations. It's exciting to reflect on the progress we're making.

Aaron Fox, the support for implicit dynamics would not have been possible without Nick!

One downside of implicit dynamics might be that it complicates the calculation of the solution's accuracy. I'm not talking about the error in the constraints reported by IPOPT; I'm talking about the ODE error. Right now, Moco does not compute or provide this error. This would be great to add in the future, and we've talked with Sam Brockie about this.

Aaron (and others): if you have suggestions for how we can improve the documentation of Moco, especially regarding how we implement implicit multibody dynamics, etc., please let us know!

Ross, I think ModOpReplaceMusclesWithDeGrooteFregly2016() is not working for you because it only replaces muscles in the model's ForceSet, not in <components>.

Ross, you can ignore the polynomial <length_approximation>; this is left over from last summer when Antoine tried including polynomial moment arms and muscle-tendon lengths; we have since removed this feature because it provided only a small benefit in convergence time.

There's a lot going on in this thread, so it is likely that Nick or I missed one of your questions. If so, please open a new forum thread.

[Brian Umberger]

This is a reply to Mritula:

Except for a few additions I made, everything in the Matlab version of example2DWalking is based on a C++ example Antoine Falisse created (located in Resources\Code\CPP). To provide targets for the tracking problem, Antoine used results from walking simulations in his recent paper:

https://journals.plos.org/plosone/artic ... ne.0217730

So, the coordinate trajectoires in referenceCoordinates.sto are synthetic data generated using the same model. That is also true for the GRFs in referenceGRF.sto. In most real applications these would be experimental data from subjects, but this approach was fine for demonstrating how to set-up and solve tracking problems.

For the predictive case, the parameters you need to specify depend on the problem you are trying to solve. For example, if the cost function is minimum-effort and you specify a fixed speed and fixed final time, then the problem is to find the minimum-effort gait for a particular speed and stride frequency. However, if you prescribe the speed but keep the final time free, then the problem is to find the optimal stride frequency that corresponds to the minimum-effort gait for a particular speed.

Thanks for your questions -- I hope this helps.

Best,
Brian

[Karthick Ganesan]

Dear Ross Miller,
Thanks for sharing your model and simulations. I tried the model with muscle wrapping in Windows 7 with Matlab 2018a. It didn't work. I am not sure of Windows 10.
Aside, How the solution N100j which is used as initial guess was obtained?
Thanks & Regards,
Karthick.

[Karthick Ganesan]

Dear Ross Miller and Aaron Fox,
This is related to mesh refinement. Betts, in his book, says that for solving sequence of optimal control problems, SQP is a better choice, mainly since barrier method alters the initial guess to lie strictly within bounds (as I recollect). Ton van den Bogert's 2011 paper also reports similar results. Could you please comment on this and/or share your experience of using mesh refinement with IPOPT?

Thanks & Regards,
Karthick.

[Ross Miller]

Hi Karthick,

By "didn't work" do you mean that it crashed Matlab, or it didn't even start running (e.g. gave some error), or did the optimization start running and iterating but just didn't converge? I've not tested the model with wrapping at all (since I'm not able to run it).

The N100j solution was produced by running that same code but having Moco generate its own initial guess (solver.createGuess) for a 20-interval problem, then I used that 20-interval solution as the guess for a 100-interval problem.

I've never used an SQP method for direct collocation problems so I can't comment on how well mesh refinement works there. With IPOPT and problems where I eventually want a solution on a "fine" grid (e.g. 100 intervals over a stride of walking) and don't have a good initial guess (e.g. a previous solution with low constraint violations), I've usually gotten good results by first solving the problem on a course grid (~10-20 intervals) then using that solution as the initial guess for the 100-interval problem. If I start right away on the fine grid with a bad guess, IPOPT convergence is spotty and it's hard to tell when/if it's going to converge.

The problem on the course grid usually either converges fast, or diverges fast and you know then to stop IPOPT and try a different initial guess or that there's a problem in the model, without waiting hours for the same bad result on the fine grid.

I've not read the Betts text recently but I think he's referring to convergence speed there. IPOPT will perturb your initial guess away from any boundaries (e.g. low muscle excitations, which we usually have) and this will mess with your initial constraint violations a lot, and it will take a while for it to get back to a solution that has both a low cost and low constraint violations. I used to think this was some mistake in my code or in IPOPT itself until Brian explained it recently =) But it still seems to handle mesh refinement well (well = converges on good solutions). SQP might do this faster.

Hope this helps,
Ross

[Ross Miller]

Also, I uploaded a 3-D version of the model and code today, if anyone is interested in using that:

https://simtk.org/projects/umocod

It seems to be working well. It still has no muscle wrapping.

Ross