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

Add mooring #107

Merged
merged 59 commits into from
Jun 20, 2023
Merged

Add mooring #107

merged 59 commits into from
Jun 20, 2023

Conversation

mabelzhang
Copy link
Collaborator

@mabelzhang mabelzhang commented Nov 22, 2022

Addresses #15

Opening for visibility. Test isn't complete yet.

Adds mooring line, chain, and seafloor.

2022-11-12-013357_2000x1690_scrot_chain

2022-11-12-015245_2000x1690_scrot_chainJoints

2022-11-12-01-30_chainOnSeafloor_scale_4x.mp4

Test WIP

The test isn't able to push the buoy all the way to the bounds to test the bounds yet. Because the heave cone is so heavy, a very large force has to be applied to move the buoy just a few meters. After some iterations, an out of bounds error happens, seeming to do with collision, and the pose values go haywire.

3000 iterations of 10,000 N (is that the unit?) in x:

[Dbg] [mooring_drift.cpp:106] iter 2275, buoy link pose: 11.9083 0.003313 -2.44625
[Dbg] [mooring_drift.cpp:106] iter 2276, buoy link pose: -2.64042e+39 -9.35589e+38 6.45126e+38
/home/developer/mbari/buoy_sim_ws/src/buoy_sim/buoy_tests/tests/mooring_drift.cpp:108: Failure
Expected: ((buoyXY - anchorXY).Length()) < (maxRadius), actual: 2.80127e+39 vs 60

ODE INTERNAL ERROR 1: assertion "aabbBound >= dMinIntExact && aabbBound < dMaxIntExact" failed in collide() [collision_space.cpp:460]
Aborted (core dumped)

Note how the pose jumps from (12, 0, -2.4) to something huge. That's not the actual test failing (looks for < 60 m from anchor), that's values going out of bounds.

There are a few things I can try:

  • Apply ocean currents (Adds support for ocean currents gazebosim/gz-sim#800) manually on the command line and see how far that pushes the buoy
  • Enable GUI in the test - is there a way to do that? I didn't find a GUI parameter in ServerConfig that we can pass to the server.
  • Write a standalone plugin (basically identical to the one in the test) to run in the simulation and watch the GUI to check things are moving correctly.

My main uncertainty is that, since I'm not seeing the GUI, even before it goes out of bounds, I don't know if the buoy is moving to the force correctly.

@mjcarroll
Copy link
Collaborator

  • Check performance against main branch
  • Check performance with collisions available
  • Test instability with ApplyLinkWrench system
  • Implement static computation to shortcut dynamics (and as something to check against).

@andermi
Copy link
Collaborator

andermi commented Dec 5, 2022

  • Implement static computation to shortcut dynamics (and as something to check against).

@mabelzhang Here's a wheel file (in a zip since github doesn't like .whl extension) for a mooring catenary calculator in python I whipped up. We can chat about it at the 2pm, or sooner if you'd like.

unzip and then pip3 install
mooring_catenary.zip

Then, you can run the gui calculator on the command line:
$ mooring_catenary_calc.py

or, import the function:
from mooring_catenary import solve_catenary
and then you can solve the catenary and compute the other values like in the calculator:

input:

V = 82.  # meters, vertical distance from buoy to anchor
H = 120.  # meters, horizontal distance from buoy to anchor
L = 160.  # meters, total length of chain
w = 20.  # N/m, weight of chain per unit length

solve:

# y: catenary curve from B to H
# B: meters, length of chain laying on the bottom, start of catenary
# c: meters, scaling factor, ratio of horizontal component of chain tension and weight of cable per unit length
y, B, c = solve_catenary(V, H, L)
T_x = c*w  # N, Horizontal component of chain tension
T_y = -w*(L - B)  # N, vertical component of chain tension @ buoy
S = L - B  # meters, length of chain suspended

output:

y = 28.14*np.cosh((x-52.38)/28.14) - 28.14
B = 52.38 (m)
c = 28.14 (m)
F_x (buoy) = -562.79 (N)
F_y (bouy) = -2152.43 (N)
F_x (anchor) = 562.79 (N)
S = 107.62 (m)

@mabelzhang mabelzhang marked this pull request as draft December 24, 2022 05:38
@mabelzhang
Copy link
Collaborator Author

mabelzhang commented Dec 24, 2022

I added the first pass of the catenary solver based on Slack discussions. Not sure if I did it correctly.
It's not working yet (hence converted PR to draft).
I'm getting solverInfo == 5 and nans, which I need to look into.
After that, the tension and force still need to be calculated from the solver's results.

@mabelzhang
Copy link
Collaborator Author

mabelzhang commented Jan 6, 2023

I went back to reread the Python diagram and the code. That made me fix a few things in the C++ solution class that I had misunderstood before.

It doesn't give me NaNs anymore, so I probably fixed something.

VSolver (in CatenarySoln.hpp. I think the fvec[0] here is supposed to be the x=H that is solved):

[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 159.431
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 161.174
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 159.431
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 161.174
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 157.793
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 165.028
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 152.096
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 189.322
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 137.247
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 4 fvec[0]: 138.85
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 159.431
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 161.174
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 5 fvec[0]: 120
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 159.431
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 161.174
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 157.793
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 165.028
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 152.096
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 189.322
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 1 fvec[0]: 137.247
[ign gazebo-1] [Dbg] [CatenarySoln.hpp:199] VSolver solverInfo: 4 fvec[0]: 138.85

This solverInfo looks okay.

HSolver (in MooringForce.cpp) just keeps printing the same thing:

[ign gazebo-1] [Dbg] [MooringForce.cpp:183] HSolver solverInfo: 5 c: 0 B: 78

This solverInfo is 5, which I think means it didn't find a solution? So maybe I have something else wrong.
Without doing anything, B is exactly L - V (160-82), which means it's finding a solution at the initial b = 0, hmm maybe that's okay? But c = 0, (since B = L - V, (L-B)^2 - V^2 is exactly 0), which I think is what we don't want? c needs to be > 0 right?

Performance is very very slow. Other than I might be doing something wrong, maybe b can be incremented faster.
I can also try moving the buoy to see if I get NaNs, or if B is still reasonable.

@andermi
Copy link
Collaborator

andermi commented Jan 6, 2023

There's an example of solver parameters to modify in ElectroHydraulicPTO.cpp

@andermi
Copy link
Collaborator

andermi commented Jan 6, 2023


// See MINPACK documentation for detail son this solver
// Parameters and defaults are (Scalar = double):
//           : factor(Scalar(100.))
//           , maxfev(1000)
//           , xtol(std::sqrt(NumTraits<Scalar>::epsilon()))
//           , nb_of_subdiagonals(-1)
//           , nb_of_superdiagonals(-1)
//           , epsfcn(Scalar(0.)) {}
  Eigen::HybridNonLinearSolver<ElectroHydraulicSoln> solver(this->dataPtr->functor);
  solver.parameters.xtol = 0.0001;
  solver.parameters.maxfev = 1000;
  solver.diag.setConstant(3, 1.);
  solver.useExternalScaling = true;       // Improves solution stability dramatically.

Try useExternalScaling=true

@mabelzhang
Copy link
Collaborator Author

With these lines locally

    vSolver.diag.setConstant(1, 1.);
    vSolver.useExternalScaling = true;

It didn't seem to make a difference to the results or the speed.

I think I must be doing something wrong in the HSolver, for it to always get solverInfo = 5... maybe I need to check my operator() again.

Tell me if this sounds wrong. The inner optimization wants y(H)-V, where y is the cosh function, so it sets fvec[0]=y(H)-V at the end. The outer optimization wants y_B(B)=H, so it first calls the inner class with initial guess x0=H, and once the inner optimization is done, it uses whatever x that is found, and just returns that? Or should it be returning y_B(B)-H (I thought that was what I was doing in f6fb01a in this line, but that was giving NaNs. Maybe I wasn't returning what I thought I was returning?

@andermi
Copy link
Collaborator

andermi commented Jan 7, 2023

but that was giving NaNs

possibly because the initial guess of B = L - V - 0 is invalid since B is strictly less than L - V

Or should it be returning y_B(B)-H

this ^^. it should be finding the B that makes x-H == 0. Just returning x tries to make x == 0 given a B which will obviously not give x == H (unless the buoy is directly over the anchor, in which case the assumptions of this catenary calc are invalid)

the assumptions are that we can use this calculation when the mooring line is pulled sufficiently taut such that the section of the line on the floor, B, is basically pulled straight with no loops/bights. Also, the catenary equation is only valid when H > B, otherwise the line just makes a straight line starting from B up to the buoy. This all means there's some radius around the anchor where we can start to use this solution.

@mabelzhang mabelzhang self-assigned this Jan 9, 2023
@mabelzhang mabelzhang mentioned this pull request Jan 10, 2023
2 tasks
@mabelzhang mabelzhang linked an issue Jan 12, 2023 that may be closed by this pull request
2 tasks
@andermi
Copy link
Collaborator

andermi commented May 13, 2023

does it still work with the old effectiveRadius code?

@mabelzhang
Copy link
Collaborator Author

mabelzhang commented May 13, 2023

Yeah 7e79a40 (right before removing effectiveRadius) works. Even 959fa9e that replaces effectiveRadius with the new toleranceL works.

The test breaks starting from the merge from main, which... there's so much new stuff. From the printouts, I can tell at least one new thing is WaveBodyInteractions.cpp, which at least wasn't printing stuff prior to the merge from main.

Maybe the test just needs to be loosened, or we can set the forces to be smaller, so that the test stops before it explodes.
Actually it fails quite quickly... it doesn't last long. It fails before any mooring printout (now it is throttled to 1 second).

If I just run the simulation, it doesn't crash. Only the test crashes.

Maybe decreasing the force in the test would help. I'd have to play with it more next week.

@mabelzhang
Copy link
Collaborator Author

mabelzhang commented May 15, 2023

1 million is too big in the force in the test.
Decreasing to 100,000 force + 2000 iterations pass the test, but the buoy is far from the max 160 m when it's stretched out.

Playing with the numbers to find a small enough number of iterations so that the buoy is pushed to near 160 without wasting too much test time:
12000 iterations at 155.006:

[Dbg] [mooring_force.cpp:116] Final buoy link pose: 235.006 99.7026 -2.21778 Final horizontal distance between buoy and anchor: 155.006

10,000 iterations at 155.119:

[Dbg] [mooring_force.cpp:116] Final buoy link pose: 235.118 99.5661 -2.74595 Final horizontal distance between buoy and anchor: 155.119

8000 iterations at 155.103:

[Dbg] [mooring_force.cpp:116] Final buoy link pose: 234.398 85.2305 -2.77727 Final horizontal distance between buoy and anchor: 155.103

5000 iterations at 153.577:

[Dbg] [mooring_force.cpp:116] Final buoy link pose: 209.852 17.9978 -1.399 Final horizontal distance between buoy and anchor: 153.577

So I think 100,000 force and 8000 iterations is a stable place to leave the test.
f68d994 does that. Will see how CI goes.

@mabelzhang
Copy link
Collaborator Author

I'm not sure why one test was missing the log, and the other test was timing out. Hard to see what's happening with a lot of printouts from WaveBodyInteractions.cpp. Maybe we should throttle that like Rhys throttled the mooring plugin printouts.

Added SDF parameter to disable plugin by default. Might as well see if the tests pass with the plugin disabled.

Signed-off-by: Mabel Zhang <[email protected]>
@mabelzhang
Copy link
Collaborator Author

This is weird. I commented out the mooring test, and the CI is still missing logs. I don't see this happening to other PRs. Have we seen this elsewhere?

@andermi
Copy link
Collaborator

andermi commented May 18, 2023

No, I haven't seen that, although github has been acting strangely the last couple days...

@andermi
Copy link
Collaborator

andermi commented May 19, 2023

I reran the test and I noticed it was running for >5hrs before it was cancelled

@mabelzhang
Copy link
Collaborator Author

Yeah I saw that. Locally, there's a lot of printouts from WaveBodyInteractions.cpp that obscures the printouts from MooringForce.cpp. I'll probably disable the printouts from the former to look into the latter more closely.

@andermi
Copy link
Collaborator

andermi commented May 20, 2023

those printouts are already masked as Debug... did you enable debug printing?

@mabelzhang
Copy link
Collaborator Author

mabelzhang commented May 20, 2023

I tried levels 1-4 in the first line of the mooring test that sets them, but running the test still prints all the debug printouts. I mean, the mooring printouts are debug too, and I want to see those.
The mooring printouts are now throttled to once per second, which is reasonable and reduces traffic, so they are getting buried by the other printouts.

@andermi
Copy link
Collaborator

andermi commented Jun 1, 2023

osrf/mbari_wec_utils#51 should fix CI test issue

@mabelzhang
Copy link
Collaborator Author

aff4053 should fix the test now, and a few bugs I found...

The plugin now has an <enable> flag that defaults to false, so the mooring is off by default.
The test calls a service to enable the mooring plugin, so that it can be tested.

That way this PR can be merged without affecting other things.
If we have time, a launch parameter to turn it on can be added in a separate PR. If not, it can be turned on by a service call to /mooring_force/enable.

@mabelzhang
Copy link
Collaborator Author

Wow, CI is finally green.

@mabelzhang
Copy link
Collaborator Author

Recap:

  • Plugin is disabled by default via the <enable>false</enable> in the model.sdf file.
  • Plugin can be enabled live by calling /mooring_force/enable service. The test does that.
  • To do in another PR: Add launch parameter to enable or disable mooring.

Copy link
Collaborator

@andermi andermi left a comment

Choose a reason for hiding this comment

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

LGTM! :)

@mabelzhang mabelzhang merged commit a4938fb into main Jun 20, 2023
@mabelzhang mabelzhang deleted the mabelzhang/mooring branch June 20, 2023 23:59
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.

Mooring
4 participants