Feature request: apply `nsubsteps` to dynamics equality constraint

[rebeccamccabe]

Feature description.
Right now nsubsteps only applies to the user-defined constraints, but I propose that it also applies to the internal dynamics equality constraint. This is needed when running with a regular wave and limiting the frequencies to a few harmonics of the excitation wave, so nfreq is rather small (ie 5), and therefore the dynamics are not enforced at very many timesteps (9). Or, allowing the user to specify nt rather than fixing it to ncomponents would work too.

If we figured out how to quantify the error as a function of number of collocation points in the pseudospectral method, perhaps in the long term nt could be calculated based on a user-specified error tolerance, which is how the step size is chosen in a regular ODE solver. But for now, any option in the API to change nt independently from nfreq would be helpful, whether it is directly or through nsubsteps.

Describe alternatives you've considered
Currently I am just setting nfreq higher than I otherwise would, in order to get enough timesteps. This is fine for now but eventually when I am doing bigger problems, the increased size of the optimization problem could interfere with performance.

Interest in leading this feature development?
This is more of an API thing so I'm not particularly interested. Perhaps it aligns with the work being done in #292.

[rebeccamccabe]

This is also probably the reason that I see weird behavior when running with nfreq=1, because the dynamics are only enforced at one timestep. The behavior is shown in the plot below, where the velocity is zero despite nonzero position and acceleration. The plot is made from the time domain results of wec.post_process and pto.post_process with nsubsteps=16. It might be helpful to have some sort of warning thrown in this case.

[rebeccamccabe]

After thinking about it some more, instead of the default nt = 2 * nfreq (or nt = 2 * ifreq_end in #292), it may make more sense to set nt = npoints * max(freqs) / min(freqs) where npoints is the desired number of points per period of the sin wave of the highest frequency (~10 minimum to get a choppy sin wave shape, ~30 to be mostly smooth) and the multiplier max(freqs) / min(freqs) accounts for how many periods of the highest frequency elapse over one cycle of the lowest frequency (which happens to equal nfreq if frequencies are evenly spaced from 0). Then users could set npoints based on desired tradeoff between accuracy and speed.

[michaelcdevin]

To recap the discussion regarding this issue today:

• Applying nsubsteps to the dynamics equations would effectively make it redundant with the frequency vector, so
• It won't be feasible to add an npoints parameter, as the number of time steps and number of frequencies are intrinsically coupled in WecOptTool due to its application of the Fourier transform in the pseudo-spectral method
• Set a different starting frequency. Closes #223 #292 will address the practicality problems of using far more frequencies than is of interest for a particular problem
• The development team will work on updating the language in the documentation and/or tutorials to make more clear the purpose of adding nsubsteps versus adding more frequencies (I will create a new issue for this)

@rebeccamccabe was there anything else you were hoping to be addressed with this issue? If not, I will go ahead and close it.

[cmichelenstrofer]

• The development team will work on updating the language in the documentation and/or tutorials to make more clear the purpose of adding nsubsteps versus adding more frequencies (I will create a new issue for this)

The main idea is this:

• The frequency array (defined by nfreq and f1) is for the dynamic equations and nsubsteps is if you want to enforce your additional constraints at more locations. This allows you to save some computation by not having to enforce the dynamics on the denser discretization. Typically you can get away with enforcing the dynamics at fewer points than the constraints.

[rebeccamccabe]

Sounds good. The main thing I’m still unclear about is the intuition on why the frequency vector and time vector are coupled in the pseudo spectral method. I looked at Giorgio’s thesis and the unnumbered equation below eq 4.96 seems to show this relationship but doesn’t explain. My intuition is saying that even with a single frequency, it’s possible to discretize the sin wave with a different number of points in time and enforce dynamics at those times, so I’m not sure why the Fourier transform demands differently.

[cmichelenstrofer]

@rebeccamccabe in that case you would have more equations (equality constraints) than unknowns and the problem is over determined. I guess you could still optimize it and you would get a best fit (based on some error measure for the objective function). But if you want to solve for an exact solution you need as many collocation points (equality constraints) as you have variables. This is not an issue for inequality constraints (e.g. max force) for which the solution is a subspace not a single point. nsubsteps should be used for inequality constraints only. We will clarify this in the documentation in #314.

[rebeccamccabe]

Ah that makes a lot of sense, thanks! Does that mean the problem is also over constrained if the user adds more non-dynamics equality constraints than there are variables in x_opt?

[cmichelenstrofer]

Yes, that is correct.