Skip to content

Attempt at solving a simple problem with friction using a slack variable. #476

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: master
Choose a base branch
from

Conversation

moorepants
Copy link
Member

@moorepants moorepants commented May 31, 2025

Attempt to fix #475

Observations:

  • The solutions for just going right or just going left work as expected but require good initial guesses.
  • It is very difficult to get the exact expected solution at the friction switching point.
  • I often get "error in step computation" from IPOPT, but don't recall seeing this in other problems in the past. Different initial guesses will sometimes cause an error. This will occur in the Cython or NumPy backends differently for the same problem conditions.
  • It seems to require a very good initial guess.

@moorepants
Copy link
Member Author

Example 100 node solution:

image

@moorepants moorepants marked this pull request as ready for review June 1, 2025 07:51
@moorepants
Copy link
Member Author

@Neville-N this may be helpful for you as an option to have elements in your equations that switch discontinuously (in alternative to making smoothened step functions).

@Peter230655 I'm curious what you can do with this and whether you can get better solutions than I am. Have a look!

m*v(t).diff(t) - Ffp(t) + Ffn(t) - F(t),
# following two lines ensure: psi >= abs(v)
psi(t) + v(t), # >= 0
psi(t) - v(t), # >= 0
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Peter and I discussed making this psi(t) - v(t)**2 instead.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or maybe it needs to be psi - sqrt(v**2).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried this and it seems to work.

@Peter230655
Copy link
Contributor

I will definitely look at it on Tuesday. Tomorrow I will be on the road and too cumbersome on ipad.

@Peter230655
Copy link
Contributor

Example 100 node solution:

image

Just running your code I do not get the nice results you show, but totally unreasonable results.
I will check some more tomorrow.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 5, 2025

Update 5-6-25:

  • With slack variables, I still cannot get the nice results shown above. It says it converges, but Ffp, Ffn have unreasonable values initially.
  • If I model the friction as Ffp = Fr * smooth_step(v(t)) and Ffn = Fr * smooth_step(-v(t)) it converges 'immediately', and final time = 0.65 sec. for any initial conditions. Fr = mu * m * g

@moorepants
Copy link
Member Author

It says it converges, but Ffp, Ffn have unreasonable values initially.

Yes, that happens with me often in this problem.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 5, 2025

I wonder whether this slack variable approach is suitable for opty/Ipopt:
Essentially

Ffp = 0 : v >=0
    = mu * m * g:  v < 0

Same for Ffn.
This is a discontinuous change, which as per my experience opty does not 'like'.
If v = psi = 0, then Ffn, Ffp may be what they want, which is not mechanically sensible?

If I bound like 0 <= Ffp, Ffn <= mu * m * g, which is correct mechanically I think, I do not get any reasonable solution.

@moorepants
Copy link
Member Author

IPOPT requires that the objective function and the constraints for the NLP are continuous. As far as I can tell, this results in a continuous formulation of the NLP. There is no necessity that the optimal control differential equations are continuous, IPOPT doesn't even know anything about the differential equations.

@moorepants
Copy link
Member Author

moorepants commented Jun 5, 2025

Also, note that this is the same formulation we used for solving the skateboard ollie, which is also solved by IPOPT.

@moorepants
Copy link
Member Author

One more thing: all unknown inputs to this system can be discontinuous with respect to time. We have many problems with such solutions in the examples.

@Peter230655
Copy link
Contributor

Also, note that this is the same formulation we used for solving the skateboard ollie, which is also solved by IPOPT.

This is what kept me wondering, why it does not easily work here.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 5, 2025

One more thing: all unknown inputs to this system can be discontinuous with respect to time. We have many problems with such solutions in the examples.

Very true, I did not think of this!
Ffp, Ffn are discontinuous w.r.t. v(t), so only indirectly w.r.t. time.

Still, it solves rapidly if I use a smooth step function, with the slack variables I am still struggling. Even if it converges the results are considerably worse ( = longer time to finish) compared to smooth step.

@Peter230655
Copy link
Contributor

As far as opty / Ipopt are concerned, would the slack variable formulation be equivalent to:

Fr = nu * m * g
Ffp = sm.Piecewise ((Fr, v(t) < 0.0), (0.0, True))
Ffn = sm.Piecewise((Fr, v(t) > 0.0,), (0.0, True))

If you think so, I can easily try this formulation.

@moorepants
Copy link
Member Author

No, that creates a discontinuous constraint function.

@Peter230655
Copy link
Contributor

Clear.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 6, 2025

I finally got these results;
I changed the integration method to midpoint.
I bound as follows: 0 <= Ffp, Ffn <= mu * m *g
It took about 4500 iterations.
The result, approx 1.2 sec is twice as large as what I get with the smooth step, which needs around 120 iterations with (just about) any initial guess.
Around 0.6 sec, when phi = 0.0, Ffp, Ffn may be what they want, so I do not take this jumping around seriously.

Could the following be possible:
With this simple problem, it is easy to see that the slack variable method does not get the optimal solution.
Maybe with the skateboard ollie this was not so easy/impossible to see?

Question: what is the advantage of the slack variable method over the smooth step method?

image

Just info: This is what the results with the sommoth step look like:
image

@moorepants
Copy link
Member Author

Question: what is the advantage of the slack variable method over the smooth step method?

In this case, its advantage is that it can find the correct solution for a system with Coulomb friction vs. the solution with an approximation of Coulomb friction. In general, the advantage is to find solutions for differential equations that have discontinuities.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 6, 2025

Clear, thanks!
Just so I understand correctly: the each eom must be continuous, but the unknown trajectories, here Ffp, Ffn, decribed by several eoms may be discontinuous?

@moorepants
Copy link
Member Author

The only stipulation is that the NLP objective and constraint functions are continuous. This says nothing about the differential algebraic equations. Opty converts the description of the optimal control of differential algebraic equations into a NLP problem via direct collocation. If the DAE's are continuous, then the NLP problem is continuous, but if the DAE's are not continuous, then you may be able to still construct a continuous NLP. That's what this example is attempting to show.

@Peter230655
Copy link
Contributor

Slowly getting it. I will play around with it more.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 7, 2025

Update 7-6-2025:

  1. Integration method ``midpoint" gives highly oscillation results, both with slack variable and with smooth step
  2. I used the result of the smooth step as initial guess for slack variable, (and psi = abs(v)). It rapidly converged to a slightly worse result (h = 0.0065 vs. h = 0.0064). At midpoint, when v changes its sign from + to - Ffn jumps to around -6 then becomes 0 as it should. Slightly perturbing this initial guess seems to move the convergence behaviour in all directions:
  • it may eliminate the negative jump of Ffn
  • it may require 14,000 iterations to reach an unreasonable result - which Ipopt still call successful.

Overall, convergence seems very sensitive to the initial guess.
Smooth step here is very insensitive to initial guesses - but as Jason pointed out, it solves a different problem.

@moorepants
Copy link
Member Author

I was wondering what it would do if given the initial guess as the smooth solution. Interesting.

Note that psi = abs(v) makes the NLP problem discontinuous. SymPy will produce a derivative for abs() but the discontinuity is present in the constraints and its Jacobian.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 7, 2025

I did not express myself correctly:
I did not change the eoms. As initial guess for psi (which does not exist in the smooth step) I used psi = v if v > 0 and psi = -v if v <0.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 7, 2025

One more thing: all unknown inputs to this system can be discontinuous with respect to time. We have many problems with such solutions in the examples.

Maybe a dumb question:
The unknown trajectory F(t) can be seleted by opty any way it wants, just limited by the bounds -400 < F(t) < 400 in this case.
F(t) is not restricted by the eoms in any way. Such unknown trajectories are generally discontinuous as selected by opty.
The unknown trajectories Ffn(t), Ffp(t) are restricted by the eoms:

Fr = mu * m * g
Ffrp(t) = Fr if v(t) < 0
        = 0 if v(t) > 0
         = any value if v(t) = 0

Similar for Ffn(t)
So are F(t) and Ffx(t) of the same 'kind' ?

@moorepants
Copy link
Member Author

In the NLP formulation there are the unknown variables x that are subject to c_l < c(x) < c_u and x_l < x < x_u. So, in that view all unknown variables are the same "kind".

@Peter230655
Copy link
Contributor

Makes sense.
I still do not understand why with slack variables it has such difficulties to converge and with smooth step it converges so easily.
The problems are different, but are they so different physically?

@moorepants
Copy link
Member Author

I think I must have something wrong in the formulation that causes the difficulty to converge.

@Peter230655
Copy link
Contributor

It is the same formulation used for the skateboard ollie, which converged easily?

@moorepants
Copy link
Member Author

It is the formulation shown in the paper, but I don't know if that is what is in the code.

@Peter230655
Copy link
Contributor

Right now I have no ideas, will continue to play around.
While I have no ideas if you could let me know about #464 I could finish it.

@moorepants
Copy link
Member Author

I chatted with Ton about this today and he pointed me to this: https://en.wikipedia.org/wiki/Linear_complementarity_problem, which points to this: https://en.wikipedia.org/wiki/Slack_variable

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 19, 2025

I think, in this sense I may have used slack variables in some of my simulations, without being aware of it:
I used a variable, say close(t) which only increased if my particle was close to a gate, like close(t).diff(t) = hump(x, a, b)
with hump(x, a, b) = 1.0 for a < x < b, zero otherwise, x the coordinate of the particle. Of course: close(t0) = 0.0

I needed this variable to be positive and tf to ensure the particle was close to the gate at some time during the simulation, but no idea of the value of close(tf) - and also not really that important, just positive.
So I introduced one more variable close1(t) and set close1(t) = aux * clsoe(t), and bounded, say, 0.1 < aux < 100
In the instance constraints I set close1(tf) = 1.0

This ensured, that the particle was close to the gate at some point, it also worked in my simulations.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 19, 2025

I think, in this sense I may have used slack variables in some of my simulations, without being aware of it: I used a variable, say close(t) which only increased if my particle was close to a gate, like close(t).diff(t) = hump(x, a, b) with hump(x, a, b) = 1.0 for a < x < b, zero otherwise, x the coordinate of the particle. Of course: close(t0) = 0.0

I needed this variable to be positive at the final time tf to ensure the particle was close to the gate at some time during the simulation, but no idea of the value of close(tf) - and also not really that important, just positive. So I introduced one more variable close1(t) and set close1(t) = aux * clsoe(t), and bounded, say, 0.1 < aux < 100 In the instance constraints I set close1(tf) = 1.0

This ensured, that the particle was close to the gate at some point in time.Iit also worked in my simulations.

@Peter230655
Copy link
Contributor

Of course I could have set close1(t) = close(t) + aux and bound -100 < aux < 0.9, if I had wanted close(tf) >= 0.1, say.
I did not think of this.

@Peter230655
Copy link
Contributor

Peter230655 commented Jun 30, 2025

I have played around with this a few hours again.
Given a good approximation, the solution of the hump differentiable approximation to the problem, it (mostly) converges, but the result is usually worse ( in the sense it takes longer to reach the final point) than hump.
Of course, it could the that the the hump approximation does not reflect the true situation well enough, but the trajectories of the slack variable solution do not look physically reasonable.

This plot shows that the 'problem' happens when v = 0 and (hence) psi = 0. I this case, Ffp, Ffn may be whatever they want.
grafik
If I bound 0 <= Ffp, Ffn <= m * mu * g it does not converge.
In summary: no progress.

@Peter230655
Copy link
Contributor

Peter230655 commented Jul 1, 2025

Update:
I managed to get a reasonable looking solution using slack variables.

  • a small change in the eoms
  • added bounds on FFp, Ffn (before it did not converge with these bounds)
  • still use the solution of the hump solution as a good initial guess
    image

Still, the convergence is quite 'sensitive' to changes in the parameters.
@moorepants : Should I push to examples-gallery, or is this still not what we want?

@moorepants
Copy link
Member Author

Nice!

I'm curious if it solves without having to give the hump solution as the initial guess. If so, I'd prefer that in the gallery to keep it minimal.

If we add this to the gallery, I'd prefer if my commits are retained and you open a PR that is based on a branch off of this PR (with your commits added).

@moorepants
Copy link
Member Author

a small change in the eoms

What is this change?

@Peter230655
Copy link
Contributor

Sorry, I did not change the eoms. I changed bounds:
I made (mu*m*g - Ffp(t) - Ffn(t))*psi(t), an equality, before it was bound to be positive.
I added bounds 0 <= FFn, FRp <= m * g * mu

@Peter230655
Copy link
Contributor

Peter230655 commented Jul 2, 2025

Nice!

I'm curious if it solves without having to give the hump solution as the initial guess. If so, I'd prefer that in the gallery to keep it minimal.

I think, it needs very good initial guesses. With your original guess, it converges to this 'solution':
image
Even with a good initial guess the solution with the slack variables is a bit worse than the hump solution (in the sense h is a bit larger), but maybe this is because hump solves a different problem.
2.
Based on my playing around, it is sensitive to many conditions, e.g. if num_nodes (called N in your original program) is large it may not converge with a good initial condition.
The hump version in theis case seems insensitive, it converges to the same solution for almost any random initial guess - but as you said, it solves a similar, but different problem.

I think, if at all we want to add it to examples gallery, we should do so with the hump as initial guess.

If we add this to the gallery, I'd prefer if my commits are retained and you open a PR that is based on a branch off of this PR (with your commits added).

Clear!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Create an example that uses slack variables for a discontinous input
2 participants