Velocity Constraints and The effect of the Scaling on the Mechanical Power

[En-Lo]

Describe the issue
We are trying to implement velocity constraints to solve for mechanical power of RM3 by following the formula:

and here is our code for the constraints:

    def const_v_pto(wec, x_wec, x_opt, waves): # Format for scipy.optimize.minimize
        v = pto.velocity(wec, x_wec, x_opt, waves, nsubsteps)
        a = x_opt[nstate_pto:]                 # auxiliary continuous time-varying decision variable, a >= 0
        return v_max - np.abs(v.flatten()) + a
    
    def const_a(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        return a
    
    def const_a_(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        return 1 - a
    
    def const_f_times_a(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        fp = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)
        return -np.abs(fp.flatten()) * a

Although the code successfully spits out the mechanical power with certain combinations of scale_x_wec, scale_x_opt, and scale_obj, these scaling factors have to be greater than 1 and influence the value of mechanical power greatly (ranging from 10^-15 to 10^-1). For instance, with scale_x_wec = 1e5, scale_x_opt = 1e2, and scale_obj = 1e3, power = -4.54e-1; with scale_x_wec = 1e1, scale_x_opt = 1e5, and scale_obj = 1e4, power = -1.22e-15.

my questions are:

1. Does our velocity constraint setup make sense?
2. Is it normal that the scaling factors affect the power output?

To Reproduce

def inner_function(i, j, k, f_max, v_max, fb, wavefreq, amplitude):

    fb.add_translation_dof(name="Heave")
    ndof = fb.nb_dofs
    
    stiffness = wot.hydrostatics.stiffness_matrix(fb).values
    mass = 208000

    f1 = 0.05 # Hz
    nfreq = 20
    freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency
    bem_data = wot.run_bem(fb, freq)
    
    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}
    
    #contraints
    nsubsteps = 4

    def const_f_pto(wec, x_wec, x_opt, waves): # Format for scipy.optimize.minimize
        f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)
        return f_max - np.abs(f.flatten())
        
    
    def const_v_pto(wec, x_wec, x_opt, waves): # Format for scipy.optimize.minimize
        v = pto.velocity(wec, x_wec, x_opt, waves, nsubsteps)
        a = x_opt[nstate_pto:]                 # auxiliary continuous time-varying decision variable, a >= 0
        return v_max - np.abs(v.flatten()) + a
    
    def const_a(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        return a
    
    def const_a_(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        return 1 - a
    
    def const_f_times_a(wec, x_wec, x_opt, waves):
        a = x_opt[nstate_pto:]
        fp = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)
        return -np.abs(fp.flatten()) * a
    
    
    ineq_cons1 = {'type': 'ineq',
                 'fun': const_f_pto,
                 }
    ineq_cons2 = {'type': 'ineq',
                 'fun': const_v_pto,
                 }
    ineq_cons3 = {'type': 'ineq',
                 'fun': const_a,
                 }
    ineq_cons4 = {'type': 'ineq',
                 'fun': const_a_,
                 }
    ineq_cons5 = {'type': 'ineq',
                 'fun': const_f_times_a,
                 }

    constraints = [ineq_cons1, ineq_cons2, ineq_cons3, ineq_cons4, ineq_cons5]
    
    
    wec = wot.WEC.from_bem(
        bem_data,
        inertia_matrix=mass,
        hydrostatic_stiffness=stiffness,
        constraints=constraints,
        friction=None,
        f_add=f_add,
        )
    
    nstate_pto = 2 * nfreq # PTO forces
    nstate_a = wec.nt * nsubsteps  #
    nstate_opt = nstate_pto + nstate_a
    
    
    #regular waves
    phase = 0
    wavedir = 0
    waves = {}
    waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)
    
    #irregular waves
    wave_cases = {
        'south_max_90': {'Hs': 0.21, 'Tp': 3.09},
        'south_max_annual': {'Hs': 0.13, 'Tp': 2.35},
        'south_max_occurrence': {'Hs': 0.07, 'Tp': 1.90},
        'south_min_10': {'Hs': 0.04, 'Tp': 1.48},
        'north_max_90': {'Hs': 0.25, 'Tp': 3.46},
        'north_max_annual': {'Hs': 0.16, 'Tp': 2.63},
        'north_max_occurrence': {'Hs': 0.09, 'Tp': 2.13},
        'north_min_10': {'Hs': 0.05, 'Tp': 1.68},
        'testing': {'Hs': 0.2, 'Tp': 5}
        }
    def irregular_wave(hs, tp):
        fp = 1/tp
        spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)
        efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, "Pierson-Moskowitz")
        return wot.waves.long_crested_wave(efth)

    for case, data in wave_cases.items():
        waves[case] = irregular_wave(data['Hs'], data['Tp'])

    obj_fun = pto.mechanical_average_power

    options = {'maxiter': 200}
    scale_x_wec = 1*10**i
    scale_x_opt = 1*10**j
    scale_obj = 1*10**k
    
    results = wec.solve(
        waves['regular'],
        obj_fun,
        nstate_opt,
        optim_options=options,
        scale_x_wec=scale_x_wec,
        scale_x_opt=scale_x_opt,
        scale_obj=scale_obj,
        x_wec_0=np.ones(wec.nstate_wec) * 1e-3,
        x_opt_0=np.ones(nstate_opt) * 1e-3,
        )

    return results.fun

length = 6
i_ = np.linspace(0, 5, length)
j_ = np.linspace(0, 5, length)
k_ = np.linspace(0, 5, length)

scale_x_wec = [i_[i] for i in range(length) for j in range(length) for k in range(length)]
scale_x_opt = [j_[j] for i in range(length) for j in range(length) for k in range(length)]
scale_obj = [k_[k] for i in range(length) for j in range(length) for k in range(length)]

t = 0
X = np.full(length**3, 0.0)
Error = np.full(length**3, np.nan)
for i in range(length):
    for j in range(length):
        for k in range(length):
            print('i, j, k', [i, j, k])
            try:
                X[t] = inner_function(i, j, k, f_max = 2000.0, v_max = 0.05, fb=RM3, wavefreq = 0.3, amplitude = 1)
            except:
                pass
            t +=1

System:

• WecOptTool version: 2.6.0
• Capytaine 2.0

[En-Lo]

Here is how we construct the RM3:

def body_from_profile(x,y,z,nphi):
    xyz = np.array([np.array([x/math.sqrt(2),y/math.sqrt(2),z]) for x,y,z in zip(x,y,z)])    # /sqrt(2) to account for the scaling
    body = cpy.FloatingBody(cpy.AxialSymmetricMesh.from_profile(xyz, nphi=nphi))
    return body

def make_RM3():
    h_f = 4.0           # height of straight section of float, before the frustum
    h_f_2 = 5.2         # height of entire float, including the frustum at the bottom
    D_s = 6.0           # diameter of spar (inner diameter of float)
    D_f = 20.0          # outer diameter of float
    mesh_density = 5

    #normal vectors have to be facing outwards
    z1 = np.linspace(-h_f_2,-h_f,mesh_density)
    x1 = np.linspace(D_s/2, D_f/2, mesh_density) 
    y1 = np.linspace(D_s/2, D_f/2,mesh_density)
    bottom_frustum = body_from_profile(x1,y1,z1,mesh_density**2)

    z2 = np.linspace(0, -h_f_2, mesh_density)
    x2 = np.full_like(z2, D_s/2)
    y2 = np.full_like(z2, D_s/2)
    inner_surface = body_from_profile(x2,y2,z2,mesh_density**2)

    z3 = np.linspace(-h_f, 0, 1+int(mesh_density/2))
    x3 = np.full_like(z3, D_f/2)
    y3 = np.full_like(z3, D_f/2)
    outer_surface = body_from_profile(x3,y3,z3,mesh_density**2)

    z4 = np.linspace(0,0,mesh_density)
    x4 = np.linspace(D_f/2, D_s/2, mesh_density)
    y4 = np.linspace(D_f/2, D_s/2, mesh_density)
    top_surface = body_from_profile(x4,y4,z4, mesh_density**2)

    RM3 = bottom_frustum + outer_surface + top_surface + inner_surface

    
    print('RM3 created')
    RM3.show_matplotlib()

    return RM3

if __name__ == '__main__':
    RM3 = make_RM3()

[dtgaebe]

Apologize the delayed response @En-Lo. This is an interesting approach that we haven't tried before. I haven't fixed all the issues and I can't confidently say that we easily can get constraints like this to work, but we can try to help you. The suggestions are implemented in the attached file
En_Issue310.zip
.

Considering the dimensions of the device, all mechanical power results are basically zero so let's ignore the scaling for now and focus on the constraints.

Issues with implementation:

1. Constraints (and dynamics) need to be formulated in the time domain
2. Amount of frequency components vs time steps (see also DOCUMENTATION: Improve explanation of nsubsteps usage #314)
3. PTO force and objective function should be adapted with non-standard optimization vector
4. Use import autograd.numpy as np instead of just numpy

Other issues:

• Numeric value for force constraint seems small for such a large device
• discontinuous jump in constraint would at least require faster dynamics (need higher frequencies, e.g through more frequencies)

Recommendations:

• Solve the unconstrained problem first
• Look at TD results and objective function
• Then introduce constraints with values that would just slightly change the solution

---

1. The auxiliary variable came directly from x_opt, which contains Fourier coefficients, so you need to transfer into the time domain, like so:

def auxiliary_nsub_td(wec, x_wec, x_opt, waves):
    aux_fd = x_opt[nstate_pto:]                 # auxiliary continuous time-varying decision variable, a >= 0
    time_matrix = wec.time_mat_nsubsteps(nsubsteps)
    aux_td = np.dot(time_matrix,aux_fd)
    return aux_td

2. I think your intention was that a(t) is a fully variable time trajectory, but you wanted more substeps than the dynamics? Thus you need the same amount of frequency components as any state (2*nfreq) and the substeps come from the time matrix as in 1.

nstate_pto = 2 * nfreq # PTO forces
nstate_a = 2*nfreq # instead of wec.nt * nsubsteps  #

3. You rightly took the last frequency components of the optimization vector, but with the "nonstandard" length of the optimization vector you should also adjust the pto force and velocity and consequently you're desired objective function.

def pto_force_td(wec, x_wec, x_opt, waves, nsubsteps=1):
    fpto_td = pto.force_on_wec(wec, x_wec, x_opt[:nstate_pto], waves, nsubsteps)
    return fpto_td
def pto_velocity_td(wec, x_wec, x_opt, waves, nsubsteps=1):
    velpto_td = pto.velocity(wec, x_wec, x_opt[:nstate_pto], waves, nsubsteps)
    return velpto_td

f_add = {'PTO': pto_force_td}

def mechpower_avg(wec, x_wec, x_opt, waves, nsubsteps=1):
    fpto_td = pto_force_td(wec, x_wec, x_opt, waves, nsubsteps)
    velpto_td = pto_velocity_td(wec, x_wec, x_opt, waves, nsubsteps)
    mechpow = fpto_td*velpto_td
    mech_energy = np.sum(mechpow)
    return  mech_energy* wec.dt/nsubsteps / wec.tf

obj_fun = mechpower_avg

the sum in the re-defined objective function requires autograd.numpy.

---

Again, interesting way to formulate a desired behavior and that you basically couple both velocity and force constraints with an entirely new optimization variable and I don't want to discourage you from finding a way to implement it, but I do wonder what the intentions are of operating a WEC like this?
Keep in mind , that WecOptTool specializes of the optimization for steady state operating conditions and while there might be no reason to limit the velocity from a generator side, I wonder if that is a desirable behavior for converting useful power?. I know, might not be you intentions, just some food for thought.
Constraining the instantaneous generator power is another option to couple force and velocity constraints. See for example e.g (23) here:
https://www.researchgate.net/publication/374138421_Incorporating_Empirical_Nonlinear_Efficiency_Into_Control_Co-Optimization_of_a_Real_World_Heaving_Point_Absorber_Using_WecOptTool

[cmichelenstrofer]

for the scaling... you should not use a single scalar for the scale_x_opt since it contains two very different quantities. Instead you should pass a vector of the same length as x_opt with the first half containing some scaling factor and the second half containing a different scaling, e.g. scale_x_opt = [a, a, a, ..., b, b, b, ...].

[jtgrasb]

@En-Lo We are working to close some old issues. Are you still working on this problem? Did @dtgaebe and @cmichelenstrofer 's suggestions help?