Skip to content
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

Mississippi steamboat on a lake #208

Merged
merged 10 commits into from
Oct 17, 2024

Conversation

Peter230655
Copy link
Contributor

@Peter230655 Peter230655 commented Aug 20, 2024

A boat is on a lake. It has a water wheel on each side to drive and to stir it, hence the name Mississippi. The control variables are the torques applied to each wheel. The goal is to get it from A to B. the objective function is a linear combination of minimizing the energy and minimizing the time needed. the drag force on a body with speed $\bar v$ is modeled as
$F_D = -c \cdot A \cdot | \bar v |^2 \cdot \dfrac{\bar v}{| \bar v |}$, where c is the drag coefficient and A is the surface facing the flow.
I do not print the EOMs as they have around 850 operations. Else, I hope all close to examples-gallery standards.

@moorepants
Copy link
Member

Error:

/home/moorepants/src/opty/docs/examples/plot_Mississippi_steamboat.rst:307: WARNING: Title underline too short.

Set up the Optimization Problem and Solve it.
--------------------------------------------
/home/moorepants/src/opty/docs/examples/plot_Mississippi_steamboat.rst:307: WARNING: Title underline too short.

Set up the Optimization Problem and Solve it.
--------------------------------------------

# - :math:`A_S`: body fixed frame of the steamboat
# - :math:`A_{LW}`: body fixed frame of the left wheel
# - :math:`A_{RW}`: body fixed frame of the right wheel
# - :math:`Dmc_S`: center of the steamboat
Copy link
Member

Choose a reason for hiding this comment

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

Is the "mc" supposed to be a subscript?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Is the "mc" supposed to be a subscript?

No, Sir. I always use, if possible Dmc to be the ceter of mass of some body. I must have copied this from some example a long time ago. If I have Dmc in a name in my code I know it is a center of mass of something.

Copy link
Member

Choose a reason for hiding this comment

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

I'm asking more because the math renders awkwardly, not because I don't like your naming conventions. You have these symbols as proper math, so I would expect them to have proper super and subscripts.

Copy link
Contributor Author

@Peter230655 Peter230655 Aug 23, 2024

Choose a reason for hiding this comment

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

I'm asking more because the math renders awkwardly, not because I don't like your naming conventions. You have these symbols as proper math, so I would expect them to have proper super and subscripts.

If I understand you correctly, something like $Dmc_S$ looks awkward?
Would DmcS be better? (I do not want to give up Dmc)

Copy link
Member

Choose a reason for hiding this comment

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

I would expect something more like:

$$ D_{{mc}_s} $$

or

$$ D^{mc}_s $$

Copy link
Member

@moorepants moorepants Aug 23, 2024

Choose a reason for hiding this comment

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

We typically use $A_o$ or $A^*$ for the mass center of rigid body $A$.

Copy link
Contributor Author

@Peter230655 Peter230655 Aug 23, 2024

Choose a reason for hiding this comment

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

We typically use A o or A ∗ for the mass center of rigid body A .

I thought capital A was 'reserved' for body fixed frames?
So, I could change to $D_{S}^{mc}$ ?

Copy link
Member

@moorepants moorepants Aug 23, 2024

Choose a reason for hiding this comment

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

You can use any letter you want for a frame, I meant that if you have a frame $X$, then the mass center is $X_o$ or $X^*$.

Copy link
Contributor Author

@Peter230655 Peter230655 Aug 23, 2024

Choose a reason for hiding this comment

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

You can use any letter you want for a frame, I meant that if you have a frame X , then the mass center is X o or X ∗ .

So, if the body fixed frame of my steamboat is $A_S$, then $A^o_S$ could be its mass center? This makes sense!
In the code AS would be the frame and AoS the center of gravity. Should I change the designation like this?

# - :math:`FP_{RW}`: point where the force due to rotation of the right
# wheel is applied

N, AS, ALW, ARW = sm.symbols('N, AS, ALW, ARW', cls = me.ReferenceFrame)
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
N, AS, ALW, ARW = sm.symbols('N, AS, ALW, ARW', cls = me.ReferenceFrame)
N, AS, ALW, ARW = sm.symbols('N, AS, ALW, ARW', cls=me.ReferenceFrame)

PEP8: no space around kwarg equals signs.

# Set up the drag forces acting on the boat.
#
# The drag force acting on a body moving in a fluid is given by
# :math:`F_D = -c \cdot A \cdot | \bar v|^2 \cdot \dfrac{\bar v}{| \bar v|}`,
Copy link
Member

Choose a reason for hiding this comment

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

Drag force is more typically:

$$ F = -\textrm{sgn}(v)\frac{1}{2} \rho C_D A v^2 $$

where $\rho$ is the density of the fluid.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

In my article, which I will remove and replace, there was the 'suggestion' to lump $\rho, C_D, \frac{1}{2}$ into one variable.
Hence I did it.

# vector of the body and :math:`A` is the cross section area of the body facing
# the flow. This may be found here:
#
# https://courses.lumenlearning.com/suny-physics/chapter/5-2-drag-forces/
Copy link
Member

Choose a reason for hiding this comment

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

You can remove this link, as it is a textbook equation. Or better link to a wikipedia article about it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You can remove this link, as it is a textbook equation. Or better link to a wikipedia article about it.

this one o.k.?
https://en.wikipedia.org/wiki/Drag_(physics)

Copy link
Member

Choose a reason for hiding this comment

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

There is a whole page about the "drag equation": https://en.wikipedia.org/wiki/Drag_equation

# ( :math:`\circ` denotes the inner product)
#
# :math:`F_{D_x} = -c \cdot A \cdot (A.x \circ \bar v)^2 \cdot \operatorname{sgn}(A.x \circ \bar v) \cdot A.x`
# :math:`F_{D_y} = -c \cdot A \cdot (A.y \circ \bar v)^2 \cdot \operatorname{sgn}(A.y \circ \bar v) \cdot A.y`
Copy link
Member

@moorepants moorepants Aug 23, 2024

Choose a reason for hiding this comment

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

I typically use \cdot for dot product and do not use a symbol for scalar multiplication, e.g. $a\hat{b}_x\cdot\hat{a}_x$.

Also, sympy prints unit vectors like: \hat{\mathbf{a}}_x.

# As an (infinitely often) differentiable approximation of the sign function,
# I will use the fairly standard approximation:
#
# :math:`\operatorname{sgn}(x) \approx \tanh( \alpha \cdot x )` with :math:`\alpha \gg 1`
Copy link
Member

Choose a reason for hiding this comment

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

If you specify all the variables to have correct assumptions (for example real=True or nonegative=True, then you can use Abs() and it will differentiate once and print with the C printers. I don't think you have to use this tangent function.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If you specify all the variables to have correct assumptions (for example real=True or nonegative=True, then you can use Abs() and it will differentiate once and print with the C printers. I don't think you have to use this tangent function.

Mathematically abs(x-a) is not differentiable at x=a, or so I understand it. How does sympy go around this?
(I had a problem with opty and abs before, but I do not think, I declared the variables to be real)

Copy link
Member

Choose a reason for hiding this comment

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

In [1]: import sympy as sm

In [2]: t = sm.symbols('t', real=True, nonnegative=True)

In [3]: v = sm.Function('v', real=True)(t)

In [4]: f = sm.Abs(v)

In [5]: f
Out[5]: Abs(v(t))

In [6]: f.diff(t)
Out[6]: sign(v(t))*Derivative(v(t), t)

In [7]: f.diff(v)
Out[7]: sign(v(t))*Derivative(v(t), v(t))

Copy link
Contributor Author

@Peter230655 Peter230655 Aug 23, 2024

Choose a reason for hiding this comment

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

Thanks!
This means, sympy defines $\dfrac{d}{dt} abs(x(t)) = 0$ at x(t) = 0, it seems.
Out[7] is not differentiable at v(t) = 0.
Would opty not 'prefer' functions whose derivatives are smooth, or does this not affect convergence?
(my gut feeling is that the more often the functions are differentiable the better the convergence - but my gut feeling was off before!!)

Should I try with sm.Abs(....)?
You wrote; t = sm.symbols('t', real=True, nonnegative=True)
I recall you told me it is better to use t = me.dynamicsymbols._t

Copy link
Member

Choose a reason for hiding this comment

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

Probably just leave it as the tanh function.

I recall you told me it is better to use t = me.dynamicsymbols._t

That is required if you use the default symbol for t in sympy.physics.mechanics. But you can set t to any variable you want.

@Peter230655
Copy link
Contributor Author

Peter230655 commented Aug 23, 2024

I (hopefully) made all the corrections you suggested, except I did not dare replacing the tanh(x) expression with the equivalent abs(x) expression. opty does not 'like' functions which are not differentiable, you told me at least once before.

NB: Let me look at my other two PRs and correct what you suggested here, before you waste your time telling me the same mistakes which I made here.

@moorepants
Copy link
Member

Here is a trick for using Abs() (but not required for this example). https://github.com/csu-hmc/gait2d/blob/master/pygait2d/segment.py#L444

@Peter230655
Copy link
Contributor Author

I used the link to drag force you gave.
I used $D^{mc}_X$ to designate the center of gravity of body X

@Peter230655
Copy link
Contributor Author

I renamed the mass centers as follows:
if $A_S$ is the body fixed frame of S, then its center of mass is $A^o_S$

@Peter230655
Copy link
Contributor Author

Here is a trick for using Abs() (but not required for this example). https://github.com/csu-hmc/gait2d/blob/master/pygait2d/segment.py#L444

Now I remember, clever! I believe you showed this no me in another context, it was this trick:
penetration = (abs(y_position) - y_position) / 2.
you showed to me then.

Would it be a sound rule of thumb for opty to avoid function which are not differentiable at certain points unless one is sure, these points will never be close to the optimum solution of opty?

@Peter230655
Copy link
Contributor Author

Hold off with this one, please.
I may have oversimplified the moment caused by a rotation of the steamboat.

@Peter230655
Copy link
Contributor Author

These two commits are the same, I just had some problem (I am still not very perfect with this github stuff), so I committed again.
All I changed was a more realistic formula for the torque ehrn the boat rotates.

I made some error in github!
@moorepants moorepants merged commit f57209e into csu-hmc:master Oct 17, 2024
17 checks passed
@moorepants
Copy link
Member

Thanks, merged!

@Peter230655
Copy link
Contributor Author

Thanks, merged!

Thank you! :-)
So, even proper names, like Mississippi should be lower case in the name.

@moorepants
Copy link
Member

I just follow PEP8 on Python module filenames:

Modules should have short, all-lowercase names. Underscores can be used in the module name if it improves readability. Python packages should also have short, all-lowercase names, although the use of underscores is discouraged.

@Peter230655
Copy link
Contributor Author

https://google.github.io/styleguide/pyguide.html#3164-guidelines-derived-from-guidos-recommendations

Clear!
So modules like this: this_is_a_module
Packages like this: this-is-a-package

@moorepants
Copy link
Member

No, they all use underscores. If you name a file a-b-c.py then you can't import it import a-b-c because the - is interpreted as subtraction in python.

@Peter230655
Copy link
Contributor Author

No, they all use underscores. If you name a file a-b-c.py then you can't import it import a-b-c because the - is interpreted as subtraction in python.

so, either thisisapackage or this_is_a_package

@moorepants
Copy link
Member

yes

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

Successfully merging this pull request may close these issues.

2 participants