BUG: Velocity and acceleration inconsistent, despite small residual

[rebeccamccabe]

Describe the bug There are some combinations of inputs (I haven't been fully systematic but it seems like it happens when there are no user defined constraints and I'm optimizing for mechanical not electrical power) that create strange high frequency oscillations in the force and acceleration signals, while the position and velocity signals do not have these oscillations. The existence of the spiky oscillations is a minor annoyance, but the bigger problem is that the velocity and acceleration are inconsistent (derivative condition not adequately enforced), even though the residual to the dynamics equality constraint is quite small.

To Reproduce

import numpy as np
import matplotlib.pyplot as plt
import sys
import xarray as xr

# switch to using local WecOptTool instead of conda package
sys.path.insert(1,'C:/Users/rgm222/Documents/Github/SEA-Lab/WecOptTool')
import wecopttool as wot

nfreq = 11
ndof = 1
nsubsteps = 1
w = np.array( [1.] )
ws = np.arange(w, w*(nfreq+1), w)
freqs = ws / (2 * np.pi)
B_h = 1/3 * w
B_hs = 1/3 * ws
K_h = np.array( [25/9 ])
m = np.array([ 1. ])
F_h = np.ones_like(ws)
name = ["PTO_Heave",]
kinematics = np.eye(ndof)    
controller = None
f_max = None

def const_f_pto(wec, x_wec, x_opt, waves):
    f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)
    return f_max - np.abs(f.flatten())
ineq_cons1 = {'type': 'ineq',
              'fun': const_f_pto,
              }
constraints = [ineq_cons1] if f_max is not None else None

nstate_opt = 2*nfreq

# analytical impedance match solution
K_p = m * w**2 - K_h
B_p = B_h
zeta = B_h / (m * w)
X_unsat = F_h[0] / (m * w**2) / (2*zeta)
Fp = np.sqrt( (B_p * w)**2 + K_p**2) * X_unsat
Fp_phase = np.arctan(B_p*ws / K_p)
real_part_Fp  = -Fp * np.cos(Fp_phase)
imag_part_Fp  = -Fp * np.sin(Fp_phase)
x_opt_0 = np.concatenate([real_part_Fp,imag_part_Fp])

loss = None
pto_imp =  np.array([[0,-1],[-1,0.01]])
pto_impedance = np.tile(np.expand_dims(pto_imp,2), [1,1,nfreq])
pto = wot.pto.PTO(ndof, kinematics, controller, pto_impedance, loss, name)

f_add = {'PTO': pto.force_on_wec}

# define impedance and reshape
impedance = m * 1j * ws + B_hs + K_h/(1j * ws)
impedance = np.reshape(impedance,(nfreq,ndof,ndof))
K_h = np.reshape(K_h,(ndof,ndof))

def make_xarrays(ws, F_h, impedance):
    # make xarrays
    freq_attr = {'long_name': 'Wave frequency', 'units': 'rad/s'}
    dir_attr = {'long_name': 'Wave direction', 'units': 'rad'}
    dof_attr = {'long_name': 'Degree of freedom'}
    dof_names = ["Pitch",]

    directions = np.atleast_1d(0.0)

    dims_exc = ('omega', 'wave_direction', 'influenced_dof')
    coords_exc = [
        (dims_exc[0], ws, freq_attr),
        (dims_exc[1], directions, dir_attr),
        (dims_exc[2], dof_names, dof_attr),
    ]
    attrs_exc = {'units': 'N/m', 'long_name': 'Excitation Coefficient'}
    exc_coeff = np.expand_dims(F_h, axis = [1,2])
    exc_coeff = xr.DataArray(exc_coeff, dims=dims_exc, coords=coords_exc,
                            attrs=attrs_exc, name='excitation coefficient')

    dims_imp = ('omega', 'radiating_dof', 'influenced_dof')
    coords_imp = [
        (dims_imp[0], ws, freq_attr), 
        (dims_imp[1], dof_names, dof_attr),
        (dims_imp[2], dof_names, dof_attr),
    ]
    attrs_imp = {'units': 'Ns/m', 'long_name': 'Intrinsic Impedance'}
    impedance = xr.DataArray(impedance, dims=dims_imp, coords=coords_imp, 
                             attrs=attrs_imp, name='Intrisnic impedance')

    return exc_coeff, impedance

exc_coeff, impedance = make_xarrays(ws, F_h, impedance)

wec = wot.WEC.from_impedance(
    freqs=freqs,
    impedance=impedance,
    exc_coeff=exc_coeff,
    hydrostatic_stiffness=K_h,
    constraints=constraints,
    f_add=f_add
)

phase = 30
wavedir = 0
amplitude = 1
waves = wot.waves.regular_wave(w/(2*np.pi), nfreq, w/(2*np.pi), amplitude, phase, wavedir)
wave = waves[:,:,0]

obj_fun = pto.mechanical_average_power

# use the analytical unsaturated solution as a guess
dc_pos = np.array([0])
real_part_pos = np.full_like(ws, X_unsat)
imag_part_pos = np.zeros_like(real_part_pos)
x_wec_0 = np.concatenate([dc_pos, real_part_pos, imag_part_pos[:-1]])

results = wec.solve(
    waves,
    obj_fun,
    nstate_opt,
    optim_options={'maxiter': 200},
    scale_x_wec=1,
    scale_x_opt=1,
    scale_obj=1,
    x_wec_0=x_wec_0,
    x_opt_0=x_opt_0,
    use_grad=False
)

x_wec, x_opt = wec.decompose_state(results[0].x)
r = wec.residual(x_wec, x_opt, wave)
print('Residual: ',r)

res_wec_fd, res_wec_td = wec.post_process(results[0],wave,nsubsteps=nsubsteps)
res_pto_fd, res_pto_td = pto.post_process(wec,results[0],wave,nsubsteps=nsubsteps)

plt.figure()
res_wec_td.pos.plot()
res_wec_td.vel.plot()
res_wec_td.acc.plot()
res_pto_td.force.plot()
res_pto_td.power.sel(type='elec').plot.line('--',x='time')
plt.legend([res_wec_td.pos.long_name, res_wec_td.vel.long_name, 
            res_wec_td.acc.long_name, res_pto_td.force.long_name,
            res_pto_td.power.long_name])

Observed behavior Run code as-is:

Optimization terminated successfully    (Exit mode 0)
            Current function value: -0.37499792016641104
            Iterations: 31
            Function evaluations: 1395
            Gradient evaluations: 31
Residual:  [ 2.06474837e-11  8.57669491e-12  2.54587462e-11  1.12005960e-10
 -5.64441827e-11 -1.41357148e-10  6.49347243e-12  6.20201668e-11
 -3.87867516e-11 -2.99227310e-11  2.62758704e-11  1.36690659e-12
 -2.11874962e-11 -1.13708154e-10 -9.00999275e-11 -8.01048117e-11
 -4.95465891e-11  1.05045750e-10 -1.82232007e-11  1.31131550e-10
  2.84490209e-11  4.59241534e-11]

The optimization converges and the residual is small, but the discrepancy between velocity and acceleration is a problem. I also plotted the three component forces, and the excitation force is smooth, while the PTO and intrinsic impedance forces have the spikes.

Also, this version is very sensitive to initial conditions, and using a different x0 results in a different (still spiky and still not obeying dynamics) solution:

Run code with obj_fun = pto.average_power: So optimizing for electrical power apparently solves the problem.

Run code with f_max = np.array( [0.94] ) So adding a force limit also apparently solves the problem.

Expected behavior When optimizing only for mechanical power with no force limit, I think that it is actually perhaps expected to have high frequency force harmonics in the PTO force, because there is nothing in the objective preventing the controller from applying them (they do not affect power). They would only go away when optimizing for electrical power, since the extra unnecessary motion would add dissipation and reduce power, or when there is a force constraint since the extra amplitude would violate the constraints. So the problem is not the presence of the spikes, but the fact that the spikes are not present in the velocity and position waveforms, indicating that the dv/dt = a portion of the dynamics is not being properly enforced, while F=ma seems to be fine.

Separately, if we wanted to get rid of the spikes, perhaps it makes sense to add a tiny norm of PTO force harmonics term to the pto.mechanical_average_power objective function, scaled such that it will have negligible effect on results and just will ensure that the additional high frequency harmonics don't occur when optimizing mechanical power. I did observe that the frequency of the spikes increases when I increase nfreq, so maybe it is only the highest harmonic where these oscillations occur.

System:

• OS: Windows 10
• Python version: 3.11.7
• WecOptTool version: main branch at 598e875cd48751ddb652657391205f6e208714f0

[cmichelenstrofer]

The residual (dynamic equations) is small, but the dynamics are wrong. I am thinking is this an issue with scaling.

Possible solution:

• change the residual output to be unscaled.

[dtgaebe]

The spikiness originates from the imaginary part of the last frequency component... I don't know why the solver converges to have content there in x_opt. When artificially removing the very last component we get what we'd expect:

x_opt[-1] = 0
x_wec[-1] = 0

res_artificial  = results[0]
res_artificial.x = np.concatenate([x_wec, x_opt])

res_wec_fd, res_wec_td = wec.post_process(res_artificial,wave,nsubsteps=nsubsteps)
res_pto_fd, res_pto_td = pto.post_process(wec,res_artificial,wave,nsubsteps=nsubsteps)

plt.figure()
res_wec_td.pos.plot()
res_wec_td.vel.plot()
res_wec_td.acc.plot()
res_pto_td.force.plot()
res_pto_td.power.sel(type='elec').plot.line('--',x='time')
plt.legend([res_wec_td.pos.long_name, res_wec_td.vel.long_name, 
            res_wec_td.acc.long_name, res_pto_td.force.long_name,
            res_pto_td.power.long_name])

[cmichelenstrofer]

might be related to #225?

[cmichelenstrofer]

This is not a bug. If you use average mechanical power as the objective function, an unstructured controller, and no constraints, then the Nyquist component of the PTO force and the WEC position is arbitrary. This is the high frequency component you are observing. If you re-run this case you should see a different magnitude for it each time. The Nyquist component of the velocity is zero which is why you don't see the high frequency oscillations.

@dtgaebe and I wrote this up showing that these components are arbitrary in this case: Power_last_frequency.pdf. We are also documenting this in #225 and closing it.

[rebeccamccabe]

The Nyquist component of the velocity is zero which is why you don't see the high frequency oscillations.

Does this mean that acceleration is not expected to match the derivative of velocity? If I want them to match, would I need to redefine the derivative matrix so that the second derivative is phase shifted 90°?

[cmichelenstrofer]

Correct, the acceleration does not match the derivative of the velocity at the Nyquist frequency. There is no way of making position-velocity-acceleration consistent at the Nyquist frequency. However if you truncate your Fourier expansion at a high enough frequency, the Nyquist component should not be significant. The problem here is that for this specific problem the Nyquist component is arbitrary and can be quiet large. The best option if you want to use average mechanical power as the objective and an unstructured controller might be to constraint the Nyquist components to be equal to zero. Here is more info on FFT-based differentiation.