Unable to capture smoothing dynamics with hydraulics style kinematics

[degoeden]

This one is a little complicated, so please let me know if there is more I can do to make this easier to follow. Not sure if I'm getting to a point where I'm going too far beyond what WecOptTool is supposed to do. I appreciate any help you are able to provide.

Describe the bug
When using hydraulics style kinematics (PTO velocity = piston area*abs(WEC velocity)) combined with smoothing dynamics (extra pole in PTO), the solver is unable to solve, unless the pole is sufficiently far away from the origin and the effects are small.

To Reproduce

1. Use the code in this gist
2. Patch the bug in the fd_to_td() function present in version 2.7.0 by following issue 321 fd_to_td() bug #329
3. Run code

Expected behavior
I expect the combination of the kinematics I introduced and the smoothing controller to find a PTO force on the WEC with this sort of shape:

Specifically I would only expect WecOptTool to show the far right side of the plot (when it settles).

Observed behavior
No convergence, final "result" has a force that looks something like this:

System:

• OS: Ubuntu 22.04
• Python version: 3.10.13
• WecOptTool version: 2.7.0 - but with the changes described in step 2 of To Reproduce
• Capytaine version: 2.0

Additional information
My guess is that the issue is occurring because the force signal on the WEC from the PTO should have a discontinuity when crossing the x axis, and this discontinuous behavior doesn't balance nicely with the sine excitation force that WecOptTool assumes the waves provide.

[cmichelenstrofer]

This is a neat case! I was able to reproduce the same results and what you did makes sense to me. The only thing I suggest double checking is that pto.force_on_wec is still doing what you want it to do (it uses the kinematics matrix). Unfortunately, it might just be a limitation of WecOptTool and Fourier bases with dealing with discontinuities.

[degoeden]

Thanks for verifying that this isn't just an error from my end. Thinking that maybe this is further motivation for #280? I will try and find another work around in the meantime.

[cmichelenstrofer]

@degoeden can you write the mathematical model you want to use? We are discussing and we don't completely follow the implementation.

[degoeden]

On the kinematics side, the math is pretty simple, but problematic. The relationship between the velocity of the WEC (v) and the velocity of the PTO, in this case hydraulic fluid flow rate (Q), is just Q = |v|*piston_area.

Similarly for the relationship between PTO force, in this case pressure (P), and the PTO force on WEC (F_pto) is just F_pto = P*-sign(v)*piston_area. The "-sign(v)" term just handles the switching of the side of the piston with pressure depending on direction of travel.

The PTO dynamics are also not terribly complex. There is simply a resistive and a reactive component. The relationship between flow rate (Q) and pressure (P) looks like this. Q = P/R + (dP/dt)*C. R represents the hydraulic resistance of the hydraulic circuit and C represents the capacitance (hydraulic accumulator). When we solve this differential eq for the transfer function we get P/Q = [C*j*omega + 1/R]^-1. Leaving us with the added pole dynamics I used in the gist, and when combined with the kinematics create what I believe to be an unsolvable problem from WecOptTool's perspective. The pressure fluctuations gets smoothed out and as a result you see discontinuities when the WEC changes direction since the pressure never drops to zero.

[dtgaebe]

Hi @degoeden ,

Here is a key things that will solve the convergence issue for the regular wave.
Get rid of the frequencies you don't need. The nfreq = 10 for Nf1 = 2 gives you 5 harmonics. Which you might need to capture the absolute value in the kinematics.

# Make a wave
amplitude = 1
wavefreq = 0.3

# set the frequency array as function of the regular wave
Nf1 = 2 #number of wavefreq in freq array (equivalent to number of period sin TD results)
f1 = wavefreq/int(Nf1)
nfreq = 10
freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency
bem_data = wot.run_bem(fb, freq)

phase = 0
wavedir = 0
waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)

I would also set the objective function scale to 1 (so that it's more intuitive)

scale_obj = 1e0

Lastly, after comparing the velocity and PTo force results, I think you should also flip the sign in the controller function.

force_fd = -H*flow_fd

I'm not sure if you can get the dynamics you expect with this model.
Where did you're plot come from?

[degoeden]

Thank you for looking into this @dtgaebe,

I've implemented your suggested edits and now also see that the solver is able to successfully terminate.

For the changes to the frequency list, I'm wondering if you could help inform me why that is a key fix. In my head more frequencies = more ability to capture non-linearities (like the absolute value function in the kinematics), so I was thinking that I may have needed even more frequencies, not less.

For the sign change, that seems to match right, I must have made an error when thinking about how the kinematics function would trace the force back to the WEC.

As for the the plot of my expectations, I've gotten the same plot from a Simscape and a Simulink model of the PTO system when I input a sine wave of velocity. I could put together a gist for those models if that would be helpful. I definitely do expect the result to look like that plot though. The kinematics of the piston combined with the accumulator (pole) should (to my understanding) create a system that effectively behaves like a full wave rectifier with a smoothing capacitor.

[dtgaebe]

Hi @degoeden ,

You're right, but you also remove a lot of frequencies that don't contain any energy content, significantly reducing the problem size. The frequencies don't contain any energy because you're using a regular wave. I hope this plot helps a bit:

The 5 harmonics I initially chose is somewhat arbitrary and stems from the experience that you can express the absolute value fairly well using a five component Fourier series.
Me wanting two periods in the results resulted in half the frequencies being irrelevant -- bad practice.
Best practice for a regular wave would actually be to just set the fundamental frequency f1 to the wave frequency.
The frequency array f1 = wavefreq, nfreq=8 would contain all the frequencies (from your initial array) that matter to the problem, for a fraction of the problem size:

For the sign change, that seems to match right, I must have made an error when thinking about how the kinematics function would trace the force back to the WEC.

Actually, I thought about your problem more, it seems like I should have flipped the signs elsewhere. I just ensured that the pto.force_on_wec opposes the wec.velocity. But since your pto.force represents the general PTO effort variable, in this case pressure, we should ensure that the PTO force in the PTO frame pto.force is positive at all times (and not negative as in the plots above)

As for the the plot of my expectations, I've gotten the same plot from a Simscape and a Simulink model of the PTO system when I input a sine wave of velocity. I could put together a gist for those models if that would be helpful. I definitely do expect the result to look like that plot though. The kinematics of the piston combined with the accumulator (pole) should (to my understanding) create a system that effectively behaves like a full wave rectifier with a smoothing capacitor.

I think the essential piece that you don't model here, but must be taken into account in your other models, is that the PTO force should only act on the WEC ones it's moving. It should only be moving ones the excitation force exceeds the opposing PTO force. This would be challenging to model in WecOptTool and is not currently supported. You might be able to find all the indexes of the collocation points where Fex<Fpto. On those collocation point you could formulate a state constraint that the velocity needs to equal zero. The tricky part here is that for every iteration you would have a different number of constraints. These dynamic constraints might require rewriting of core.py and I'm unsure of the solvers we use would allow for a dynamic constraint structure in the first place....

On a side note, I'm unsure how you can input a sine wave velocity (if it's supposed to represent a WEC), because the opposing PTO force would change the WEC's velocity, feedback....). Do you have experimental data and a velocity profile of an actuator? I think then you should be looking at the input force profile, no?