Using Moco to calculate the kinematics of a DOF/coordinate

[Sietse Achterop]

Hello list,

I understood that it should be possible to use Moco to solve problems where not all DOF's are given in that it will calculate the unknowns DOF's.
In my use case I have a model consisting of a "person" sitting on a box ("boat") that tries to push this box over the ground ("water").
The motion of the persion is described but the motion of the box should be determined.

In a simplified version the person is only a simple "leg" attached to this "boat".
The foot of the leg also is a box. The leg then can push the system over the ground given that both boxes have friction wrt to the ground.
The foot has a large and the boat has a small friction coefficient.

The pusher

pusher.png (69.52 KiB) Viewed 760 times

The leg pushes the box, then lifts its foot and moves to the position to do another push.
The idea is to prescribe the motion of the leg and let Moco determine the movement of the box.
I created an example that can be found here:

https://github.com/SietseAchterop/Rowin ... ter/Pusher

The model, Pusher_noC.osim, is calculated from pusher_nocontrollers.py.
The motion, trajectory.trc, is calculated from trajectory.py.
pusher_CMC_video.mp4 is a video that shows the pusher in action, created with CMC. See also PS below.

In calculating the motion I assumed that the boat is not moving. When using the IK-tool of Opensim a trajectory.mot is created.
The idea is to remove the coordinates that describe the position of the boat and give the rest to Moco.

I tried MocoInverse, but that does not work because it needs all coordinates.
But MocoTrack seems to need the trc-file. But that is describing the movement wrt to ground, so it will not let the boat move at all.

Or do I have to create a MocoStudy for this, and how to go about this?

Thanks in advance,
Sietse



PS. I have this "pusher" working with CMC, but more complicated versions do not work. That seems to be a limitation of CMC.
In the more complicated model the person is using its arms to handle oars, as in a rowing boat.
Hopefully, when the above simple model is working with Moco I can try more complicated models.

[Nicholas Bianco]

Hi (again) Sietse,

Rather than providing marker trajectories via a TRC file, you can specify the trajectories of some of the coordinates in the model. To do this, use setStatesReference() rather than setMarkersReference(). You can then use set_states_weight_set() to set the weights for individual tracked coordinates. Just set the boat coordinate weights to zero and you should be good to go.

Let me know if you have any further questions.

Best,
Nick

[Sietse Achterop]

Hi Nick,

I created a script pusherMocoTrack.py starting from the function muscleDrivenStateTracking from the example3DWalking/exampleMonoTrack.py. See github, it is also attached below.
I changed the model, and did not load an external load. Also I tried to set the proper weights for the coordinates.
And finally I wrote the result to mocotrack_solution.sto.
I removed the boat-coordinates from the original trajectory.mot, it now only contains the 3 angles in the leg.
I assume that the default values for the boat coordinates will be those from the model.

Execution yields a result, but the trajectory is completely wrong!

The file muscle_driven_state_tracking_tracked_status.sto that is created is almost a copy of trajectory.mot, apart from the change to radians. So far so good.
But the result stored in mocotrack_solution.sto and seen in OpenSim is completely different!

I tried to set the weights for the coordinates as follows:

Code: Select all

    coordWeights = osim.MocoWeightSet()
    coordWeights.cloneAndAppend(osim.MocoWeight("bJoint_3", 0))
    coordWeights.cloneAndAppend(osim.MocoWeight("bJoint_4", 0))
    coordWeights.cloneAndAppend(osim.MocoWeight("baseangle", 20))
    coordWeights.cloneAndAppend(osim.MocoWeight("kneeangle", 20))
    coordWeights.cloneAndAppend(osim.MocoWeight("bladeangle", 20))
    track.set_states_weight_set(coordWeights)

I did not quite know whether I needed the customization as in the example, and tried to set the WeightForControl of the boat-actuators to 0, but nothing changed.
There is no feedback in the std-output about the (non-default) setting of the coordinate- or force- weights. Is that OK?

Isn't a bit strange that the solver claims success where the trajectory is so different from the goal?

Thanks again, Sietse

PS. I tested on a 4 and a 12 core machine. The first needed 266 iterations, and the second about 61 for the same program. Shouldn't that be about the same?
Or is it counting something different?

[Nicholas Bianco]

Hi Sietse,

The solver "success" status means that the optimization converged, but the converged solution may not necessarily be the one you are looking for. This could be due to a number of different factors, including the choice of initial guess, cost weights, solver settings, etc.

My first guess (without seeing the full solution) is to try increasing the global state tracking weight. I'm not sure what the current version of the model is, but I would start by adding torque actuators to all degrees-of-freedom and temporarily removing the custom blade forces. Get that to work, then slowly add complexity back into your model.

I would set some non-zero weight for all actuators, as this improves optimization convergence.

Regarding the number of iterations on different machines: IPOPT sometimes behaves differently on different machines (e.g., Mac versus Windows, or two different Windows versions). A small difference in the initial guess or initial search direction could have noticeable effects on optimization performance. But as long as the final solutions aren't significantly different, you should be okay.

Best,
Nick

[Sietse Achterop]

Still no luck regrettably.
I increased the global_state_tracking_weight to 100.
Also did modelProcessor.append(osim.ModOpAddReserves(250)) to add actuators, but is this actually relevant because the model has no muscles and already coordinate actuators?
Finally set all weights for the actuators to non-zero.
I did this in a number of steps, but nothing changed, also the number of iterations stayed exactly the same (376).
The current version of pusherMonoTrack.py is:

Code: Select all

import opensim as osim

def muscleDrivenStateTracking():
    track = osim.MocoTrack()
    track.setName("muscle_driven_state_tracking")
    modelProcessor = osim.ModelProcessor("Pusher.osim")
    modelProcessor.append(osim.ModOpAddReserves(250))
    track.setModel(modelProcessor)
    track.setStatesReference(osim.TableProcessor("trajectory.mot"))
    track.set_states_global_tracking_weight(100)

    coordWeights = osim.MocoWeightSet()
    coordWeights.cloneAndAppend(osim.MocoWeight("bJoint_3", 10))
    coordWeights.cloneAndAppend(osim.MocoWeight("bJoint_4", 10))
    coordWeights.cloneAndAppend(osim.MocoWeight("baseangle", 20))
    coordWeights.cloneAndAppend(osim.MocoWeight("kneeangle", 20))
    coordWeights.cloneAndAppend(osim.MocoWeight("bladeangle", 20))
    track.set_states_weight_set(coordWeights)
    track.set_allow_unused_references(True)
    track.set_track_reference_position_derivatives(True)
    track.set_initial_time(0.0)
    track.set_final_time(5.0)
    track.set_mesh_interval(0.08)
    sol = track.solve()
    sol.write("mocotrack_solution.sto")
    return

muscleDrivenStateTracking()

I updated github with this, and also the model (pusher.py, Pusher.osim) and the result mocotrack_solution.sto.

UPDATE: In this model I do not use a custom force, only ElasticFoundationForce.

Thanks again, Sietse

PS. The trajectory has a repetition rate of 30 strokes/minute. I also tried a very slow rate of 3 strokes per minute but no difference.

[Nicholas Bianco]

Hi Sieste,

The number of iterations and solution are exactly the same? If so, I'm wondering if the the MocoProblem being solved somehow isn't changing with the changes you're making in the code. You can check the summary of the MocoProblem that prints to the console to ensure that you are solving the correct problem.

Nick

[Sietse Achterop]

Hi Nick,

The trajectories at least look very much the same, and very different from the intended one.
But the console output is, apart from the timing exactly the same. Here a complete diff of 2 sessions

Code: Select all

4c4
< [info] Mon Feb 14 20:52:11 2022
---
> [info] Mon Feb 14 20:53:26 2022
7c7
<   state_tracking. MocoStateTrackingGoal, enabled: true, mode: cost, weight: 100.0
---
>   state_tracking. MocoStateTrackingGoal, enabled: true, mode: cost, weight: 10.0
507,508c507,508
< Total CPU secs in IPOPT (w/o function evaluations)   =      8.874
< Total CPU secs in NLP function evaluations           =    253.065
---
> Total CPU secs in IPOPT (w/o function evaluations)   =      8.903
> Total CPU secs in NLP function evaluations           =    258.624
512,517c512,517
< callback_fun  |  11.93ms ( 57.36us)  11.92ms ( 57.31us)       208
<        nlp_f  |  11.78 s (  8.55ms)   2.28 s (  1.65ms)      1377
<        nlp_g  |  23.31 s ( 16.93ms)   3.51 s (  2.55ms)      1377
<   nlp_grad_f  |  29.98 s (140.75ms)   3.62 s ( 16.98ms)       213
<    nlp_jac_g  | 197.68 s (517.48ms)  25.99 s ( 68.02ms)       382
<        total  | 271.00 s (271.00 s)  43.65 s ( 43.65 s)         1
---
> callback_fun  |  11.66ms ( 56.08us)  11.63ms ( 55.91us)       208
>        nlp_f  |  12.01 s (  8.72ms)   2.58 s (  1.88ms)      1377
>        nlp_g  |  23.67 s ( 17.19ms)   3.81 s (  2.77ms)      1377
>   nlp_grad_f  |  30.17 s (141.65ms)   3.61 s ( 16.93ms)       213
>    nlp_jac_g  | 202.63 s (530.44ms)  26.24 s ( 68.69ms)       382
>        total  | 276.79 s (276.79 s)  44.56 s ( 44.56 s)         1
520,521c520,521
< [info] Elapsed real time: 44 second(s).
< [info] Mon Feb 14 20:52:55 2022
---
> [info] Elapsed real time: 45 second(s).
> [info] Mon Feb 14 20:54:11 2022

In the beginning you see the difference of the weight, 10 or 100.
And in the end, in the summery only the times are a little different. For the rest, the entire 40KB files are the same.

I am running on Linux, and opensim is compiled from source. No customization.
I then tried on a Windows machine (binary install), but the weird thing is that there the solver fails after 3000 iterations.
It ends with:

Code: Select all

2999r 2.3322159e-02 7.39e+01 1.42e+02  -3.5 5.72e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
3000r 2.3459691e-02 1.67e+00 7.63e+00  -3.6 6.01e-01    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 3000

                                   (scaled)                 (unscaled)
Objective...............:   2.3459690898518099e-02    2.3459690898518099e-02
Dual infeasibility......:   7.6268774991044328e+00    7.6268774991044328e+00
Constraint violation....:   4.8764258177433750e-01    1.6735707153933941e+00
Complementarity.........:   2.3563661191793763e-04    2.3563661191793763e-04
Overall NLP error.......:   7.6268774991044328e+00    7.6268774991044328e+00


Number of objective function evaluations             = 13062
Number of objective gradient evaluations             = 2043
Number of equality constraint evaluations            = 13087
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 3041
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =     93.906
Total CPU secs in NLP function evaluations           =    248.739

EXIT: Maximum Number of Iterations Exceeded.
         nlp  :   t_proc      (avg)   t_wall      (avg)    n_eval
callback_fun  | 415.00ms (136.56us) 403.32ms (132.71us)      3039
       nlp_f  |  21.11 s (  1.62ms)  21.18 s (  1.62ms)     13062
       nlp_g  |  30.95 s (  2.36ms)  30.97 s (  2.37ms)     13087
  nlp_grad_f  |  25.21 s ( 12.33ms)  25.23 s ( 12.34ms)      2044
   nlp_jac_g  | 171.21 s ( 56.28ms) 171.18 s ( 56.27ms)      3042
       total  | 342.66 s (342.66 s) 342.66 s (342.66 s)         1

I would guess that the Windows version is the most tested one. On the other hand it looks a bit strange to me that the solver could not solve such a not too complicated model. But that's just my intuition.

I also ran it on 2 other Linux machines. One behaves like my main Linux machine and the other behaves like Windows with a failing solver.
Weird...

What can I do now?
Regards, Sietse


PS. On the windows machine I get the following warning:

Code: Select all

python.exe : Warning: intermediate_callback is disfunctional in your installation. You will only be able to use stats(). See
https://github.com/casadi/casadi/wiki/enableIpoptCallback to enable it.

But I guess that this is probably not relevant.

[Nicholas Bianco]

Hi Sietse,

After taking a look through some of the files and your current problem setup, here's a few recommendations:

1) When setting weights for the tracked coordinate values, you need to provide a full path to the coordinate values. For example,

Code: Select all

coordWeights = osim.MocoWeightSet()
coordWeights.cloneAndAppend(osim.MocoWeight("/jointset/baseJoint/baseangle/value", 20))

You can set individual weights for coordinate values and speeds:

Code: Select all

coordWeights = osim.MocoWeightSet()
coordWeights.cloneAndAppend(osim.MocoWeight("/jointset/baseJoint/baseangle/value", 20))
coordWeights.cloneAndAppend(osim.MocoWeight("/jointset/baseJoint/baseangle/speed", 1))

If you open the Moco solution files, we print a breakdown of the objective values in the file headers. In the "mocotrack_solution.sto" file, the state tracking goal was exactly zero, meaning that you weren't actually tracking anything, which explains the problems you're having.

2) The trajectories for coordinates "bJoint_3" and "bJoint_4" are all zeros, so make sure the tracking weights for these coordinates (both values and speeds) are zero. For the initial coordinate tracking problem it won't matter, but once you add back in all the contact and blade forces, you don't want to track these coordinates, as the problem will then try to force the boat to stay still.

3) For now, remove all contact forces including the ElasticFoundationForces still in the model. Get a problem to solve with only coordinate actuators first. Then, start adding complexity back to the model.

4) Try using a denser mesh, perhaps a mesh interval between 0.02 and 0.05. For example,

Code: Select all

    
track.set_mesh_interval(0.025)

5) While you're still debugging your optimization, solve the problem over a smaller time window (i.e., [0, 1] seconds rather than [0, 5] as it is now).

Let's try all of these steps and see how it works. Once we get something very simple working, we can build up to the problem that you're trying to solve.

Best,
Nick

[Sietse Achterop]

Hi Nick, and thanks again for all the effort!

With state tracking goal being zero, did you mean objective_state_tracking=0.000000?
I have never seen anything else than zero in the mocotrack_solution.sto file.

I removed the ElasticFoundationForces.
I corrected the value names, but that didn't help. Also setting values to 0, 20 or 1000 did not matter.
The mesh_interval is now at 0.025, no difference.

The trajectory still is completely wrong. Only bJoint_4 (Y-axis) is moving, the other joints don't change.
When I set gravity to 0, then bJoint_4 also stops moving and the model is completely still.
If I set the time window to say 10 seconds, then there is some movement in the leg.
Hopefully this helps.

Github is updated, and the model now is p.py -> Pusher.osim

Best,
Sietse

PS. Is there a way to give the mocotrack_solution a name? Now the first line in the .sto file seems to be missing.

[Nicholas Bianco]

Hi Sietse,

Yes, "objective_state_tracking=0.000000" either means you're tracking the data perfectly (which isn't the case) or the tracking term is all zeros.

I noticed one more thing that should fix the issue: the columns in tracking reference data also need to have the full paths to coordinate values. You can use a TableOperator to do this for you:

Code: Select all

tableProcessor = osim.TableProcessor("trajectory.mot")
tableProcessor.append(osim.TabOpUseAbsoluteStateNames())
track.setStatesReference(tableProcessor)

Hopefully this fixes it!