Set a different starting frequency. Closes #223

[cmichelenstrofer]

Closes #223

[coveralls]

Pull Request Test Coverage Report for Build 6961134641

• 0 of 0 changed or added relevant lines in 0 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage remained the same at 93.93%

---

Totals
Change from base Build 6910014540:	0.0%
Covered Lines:	2476
Relevant Lines:	2636

---

💛 - Coveralls

[cmichelenstrofer]

Status:

• The modified code works when no frequencies are skipped but converges to random (initial guess dependent) wrong solutions when even one frequency is skipped.
• We have done a thorough check and every function and gradient is as expected.

As part of debugging we tried the following in the main code:

• A: add an additional equality constraint that requires some frequency components to be zero (components that are zero anyways in the optimal solution). This showed the same behavior as the PR code... convergence to random values
• B: directly change the relevant frequency components to zero in the code for both x_wec and x_opt. Same bad behavior.
  • Going more in depth, this only breaks when the equality constraint (residual of the dynamic equations) has the frequency components in both x_wec and x_opt set to zero.

[michaelcdevin]

A: add an additional equality constraint that requires some frequency components to be zero (components that are zero anyways in the optimal solution). This showed the same behavior as the PR code... convergence to random values

@cmichelenstrofer I expanded some on this using the code below:

Python code

import autograd.numpy as np
import capytaine as cpy

import wecopttool as wot


# Frequencies
wavefreq = 0.6 # Hz
f1 = wavefreq/2
nfreq = 10
freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency

# Waves
amplitude = 0.0625 # m
phase = 30 # degrees
wavedir = 0 # degrees
waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)
extra_initial_freqs = 1
extra_initial_comps = extra_initial_freqs * 2
print(f'Wave elevations for each frequency: {np.abs(waves.loc[:,0,0])}')

# Geometry (WaveBot)
wb = wot.geom.WaveBot()  # use standard dimensions
mesh_size_factor = 0.5 # 1.0 for default, smaller to refine mesh
mesh = wb.mesh(mesh_size_factor)
fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot")
fb.add_translation_dof(name="Heave")
ndof = fb.nb_dofs

# BEM
bem_data = wot.run_bem(fb, freq)

# PTO
name = ["PTO_Heave",]
kinematics = np.eye(ndof)
controller = None
loss = None
pto_impedance = None
pto = wot.pto.PTO(ndof, kinematics, controller, pto_impedance, loss, name)
f_add = {'PTO': pto.force_on_wec}
f_max = 750.0
nsubsteps = 4

# Constraints
tol = 1e-8
def const_zero_xw_ineq(wec, x_wec, x_opt, waves):
    return tol - np.abs(x_wec[1:extra_initial_comps+1])

def const_zero_xo_ineq(wec, x_wec, x_opt, waves):
    return tol - np.abs(x_opt[1:extra_initial_comps+1])

def const_zero_xw_eq(wec, x_wec, x_opt, waves):
    return x_wec[1:extra_initial_comps+1]

def const_zero_xo_eq(wec, x_wec, x_opt, waves):
    return x_opt[1:extra_initial_comps+1]

constraints_firstfreqzero_ineq = [
                             {'type': 'ineq', 'fun': const_zero_xw_ineq},
                             {'type': 'ineq', 'fun': const_zero_xo_ineq},
                             ]

constraints_firstfreqzero_eq = [
                             {'type': 'eq', 'fun': const_zero_xw_eq},
                             {'type': 'eq', 'fun': const_zero_xo_eq},
                             ]

# WEC object
wec = wot.WEC.from_bem(
    bem_data,
    constraints=None,
    f_add=f_add,
)

wec_firstfreqzero_ineq = wot.WEC.from_bem(
    bem_data,
    constraints=constraints_firstfreqzero_ineq,
    f_add=f_add,
)

wec_firstfreqzero_eq = wot.WEC.from_bem(
    bem_data,
    constraints=constraints_firstfreqzero_eq,
    f_add=f_add,
)

# Objective function
obj_fun = pto.mechanical_average_power
nstate_opt = 2*nfreq
x_wec_0 = np.random.randn(wec.nstate_wec)
x_opt_0 = np.random.randn(nstate_opt)
options = {'maxiter': 200}
scale_x_wec = 1e1
scale_x_opt = 1e-3
scale_obj = 1e-2

# Reference case
results = wec.solve(
    waves, 
    obj_fun, 
    nstate_opt,
    optim_options=options, 
    x_wec_0=x_wec_0,
    x_opt_0=x_opt_0,
    scale_x_wec=scale_x_wec,
    scale_x_opt=scale_x_opt,
    scale_obj=scale_obj,
    )

print('============================\nREFERENCE:\n============================')
print(f'Optimal average mechanical power: {results[0].fun} W')
print(f'Optimal state vector: {results[0].x}')

# Inequality constraint case
results_ffz_ineq = wec_firstfreqzero_ineq.solve(
    waves, 
    obj_fun, 
    nstate_opt,
    optim_options=options,
    x_wec_0=x_wec_0,
    x_opt_0=x_opt_0,
    scale_x_wec=scale_x_wec,
    scale_x_opt=scale_x_opt,
    scale_obj=scale_obj,
    )

print('============================\nWITH INEQUALITY CONSTRAINTS:\n============================')
print(f'Optimal average mechanical power: {results_ffz_ineq[0].fun} W')
print(f'Optimal state vector: {results_ffz_ineq[0].x}')

# Equality constraint case
results_ffz_eq = wec_firstfreqzero_eq.solve(
    waves, 
    obj_fun, 
    nstate_opt,
    optim_options=options,
    x_wec_0=x_wec_0,
    x_opt_0=x_opt_0,
    scale_x_wec=scale_x_wec,
    scale_x_opt=scale_x_opt,
    scale_obj=scale_obj,
    )

opt_power_ffz_eq = results_ffz_eq[0].fun
print('============================\nWITH EQUALITY CONSTRAINTS:\n============================')
print(f'Optimal average mechanical power: {opt_power_ffz_eq} W')
print(f'Optimal state vector: {results_ffz_eq[0].x}')

Relevant output

Wave elevations for each frequency: <xarray.DataArray 'wave_elev' (omega: 10)>
array([0.    , 0.0625, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ,
       0.    , 0.    ])
Coordinates:
  * omega           (omega) float64 1.885 3.77 5.655 7.54 ... 15.08 16.96 18.85
    freq            (omega) float64 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0
    wave_direction  float64 0.0
    realization     int32 0
Attributes:
    units:                m
    long_name:            Wave elevation
    Wave type:            Regular
    Frequency (Hz):       0.6
    Amplitude (m):        0.0625
    Phase (degrees):      30
    Direction (degrees):  0
Optimization terminated successfully    (Exit mode 0)
            Current function value: -0.17374424555902082
            Iterations: 84
            Function evaluations: 84
            Gradient evaluations: 84
============================
REFERENCE:
============================
Optimal average mechanical power: -17.37442455590208 W
Optimal state vector: [-5.83135587e-03 -2.01622688e-06  1.64814845e-06  3.90771495e-02
 -7.31322604e-03 -1.91736455e-07  6.20508181e-07 -3.66666032e-07
  1.58528911e-06  3.47237235e-06 -9.65765026e-06 -1.82382984e-03
 -1.45634505e-03 -5.11003352e-04  2.05942038e-04  9.34183391e-06
 -9.67135827e-05  1.87657519e-02  9.86448906e-02 -1.87431786e-04
 -1.42244417e+02 -3.69525824e-02  2.38348924e-02  1.84370310e+01
 -2.41886964e+02  2.91500841e-03 -1.41675464e-02  2.14622646e-02
 -1.06618653e-01 -4.21637155e-01  1.18793723e+00  3.48520252e+02
  2.78262296e+02  1.38798060e+02 -5.60002340e+01 -3.39463004e+00
  3.52030942e+01 -8.79260860e+03 -4.62184389e+04  1.09715509e+02]
Optimization terminated successfully    (Exit mode 0)
            Current function value: -0.17374286022314755
            Iterations: 76
            Function evaluations: 76
            Gradient evaluations: 76
============================
WITH INEQUALITY CONSTRAINTS:
============================
Optimal average mechanical power: -17.374286022314756 W
Optimal state vector: [-5.83135584e-03  6.54915049e-13  5.22180600e-13  3.90770514e-02
 -7.31302392e-03  3.17679374e-07  2.01110847e-07 -6.97230712e-07
  1.56616962e-06  4.19504190e-06 -1.00783914e-05 -1.83172006e-03
 -1.46054654e-03 -5.15616200e-04  1.99121258e-04  6.47803543e-06
 -9.23157161e-05  1.87723445e-02  9.86803049e-02 -1.87431786e-04
 -1.42244416e+02  1.00000006e-08  9.99999986e-09  1.84356849e+01
 -2.41887224e+02 -7.47484074e-03 -3.77446408e-03  4.35846101e-02
 -1.05974198e-01 -5.10186923e-01  1.23997757e+00  3.50027990e+02
  2.79065039e+02  1.40051910e+02 -5.41477787e+01 -2.35250125e+00
  3.36021570e+01 -8.79569753e+03 -4.62350317e+04  1.09715509e+02]
Optimization terminated successfully    (Exit mode 0)
            Current function value: 433.0853832991693
            Iterations: 4
            Function evaluations: 21
            Gradient evaluations: 4
============================
WITH EQUALITY CONSTRAINTS:
============================
Optimal average mechanical power: 43308.538329916926 W
Optimal state vector: [-2.25606010e-01  0.00000000e+00  2.10449248e-22 -3.80185722e-02
  1.12312320e+00 -9.81625598e-01 -1.74888799e+00  3.59629816e-01
 -5.95291942e-01 -1.72934627e-01  6.16247310e-01  6.95920502e-01
  3.51975480e-01 -6.64044795e-02 -2.72168386e-01 -1.35242966e-01
  2.26920656e-01  3.62740802e-02  1.89212792e-01 -2.32769242e-03
 -5.50321331e+03 -5.83605713e-15 -9.21592893e-14 -6.76719220e+03
  1.08410166e+03  2.55233256e+04  3.66556758e+04 -2.28877868e+04
  4.04613444e+04  2.09308199e+04 -7.57769268e+04 -1.32982863e+05
 -6.72486444e+04  1.80681646e+04  7.39305764e+04  4.92134321e+04
 -8.26042516e+04 -1.69960538e+04 -8.86525368e+04  1.36254349e+03]

• Using equality constraints to set the "extra" frequency to zero leads to convergence to an incorrect value, though this value is consistent if using the same initial guesses
• Using inequality constraints to set the extra frequency to within machine precision of zero leads to the correct value
• Interestingly, adding in another constraint (e.g. the PTO force constraint from Tutorial 1) results in a "Positive directional derivative for linesearch" exception when trying to solve the equality constraint version. The other versions still solve correctly.