Using WEC Opt Tool for arrays

[NicolasBarbarin]

Hello,

I am trying to use WEC Opt Tool for arrays, for instance 3 WECs moving in heave only. That means ndof = 3. I have create the Floating Body using Capytaine and run the bem solver that works well.

Then, to create the PTO object, I am using this :

name = ["PTO_Heave","PTO_Heave","PTO_Heave"]
kinematics = np.eye(3)
controller = None
loss = None
pto_impedance = None
pto = wot.pto.PTO(1,kinematics, controller, pto_impedance, loss, name)

But I would like to have ndof = 3 and not ndof =1. When I change it in the last line of the code an error occurs... Maybe this is because of controller which has to be a 3x3 matrix. But what is the equivalent of "controller = None" for a 3 dofs-PTO ?

Then, a second issue appears. I want to define a constraint for each WEC, which mean a constraint for the 3 different degrees of freedom (heave1, heave2, heave3). How can I separate this when defining the constraint like that :

def const_pos_pto_heave(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps_p)
    return 3 - np.abs(pos.flatten())

For example I want to constraint the position of the WEC in heave to 3m. But let say I want now to constraint the postion of WEC1 to 3m and the position of WEC2 to 2m. How can I do this. I tried the following that isn't working :

def const_pos_pto_heave1(wec, x_wec, x_opt, waves):
    n = len(x_wec)//3
    x_wec_h = np.concatenate((x_wec[:n], [0]*n,[0]*n))
    pos = pto.position(wec, x_wec_h, x_opt, waves, nsubsteps_p)
    return pos_max_heave - np.abs(pos.flatten())

def const_pos_pto_heave2(wec, x_wec, x_opt, waves):
    n = len(x_wec)//3
    x_wec_h = np.concatenate(([0] * n , x_wec[n:2*n] , [0] * n))
    pos = pto.position(wec, x_wec_h, x_opt, waves, nsubsteps_p)
    return 2 - np.abs(pos.flatten())

constraints = [{'type': 'ineq','fun': const_pos_pto_heave1}, {'type': 'ineq','fun': const_pos_pto_heave2}]

This is not working because of an error of dimension.

Could you help me on this ?

Thanks in advance Nicolas Barbarin, nicolas.barbarin@corpowerocean.com

[cmichelenstrofer]

Hi Nicolas,

This should definitely be doable!

But I would like to have ndof=3 and not ndof=1. When I change it in the last line of the code an error occurs...

What error do you get? I ran that block of code and got no error (it created the PTO object no problem with ndof=3)

The separate constraints... you shouldn't need to hack it by setting part of x_wec to zero. This line:

pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps_p)

gives you the position of all three DOFs each in a different row (matrix 3xN). You have different options here. You could have three different constraints, and select the relevant row in each. Or you can have a single constraint but subtract the correct value for each element.

what is the equivalent of "controller = None" for a 3 dofs-PTO ?

controller=None use a "unstructured controller" (i.e., the control trajectory will be the numerical optimal controller). There is some discussion of this concept on the second page here: https://doi.org/10.1109/TSTE.2023.3272868.

[NicolasBarbarin]

Hi,

Thank you for your answer ! Yes, the PTO object is created correctly with no error, but then I get an error at the “wec.solve” step. Here are some more details :

My code :

[Carlos: Moved and formatted, see comment below]

Every block is running without error except the latter that returns the following :

[Carlos: Moved and formatted, see comment below]

So apparently, the error could come from the bad definition of the controller in the PTO object. But even if I try to set the controller as a 3x3 matrix as I see in the WecOptTool documentation, (firstly filled with “None” or even filled with random integers), the same error occurs…

Do you have an idea of where the problem is coming from ?

Thank you very much in advance,

Nicolas

[Carlos: Removed email reply chain]

[cmichelenstrofer]

Code from email above:

import autograd.numpy as np

#Bem_solver
f1 = 0.02
nfreq = 50
freq = wot.frequency(f1, nfreq, False)
bem_data = wot.run_bem(FB, freq)

#PTO
name = ["PTO_Heave1","PTO_Heave2","PTO_Heave3"]
kinematics = np.eye(3)
controller = None
loss = None
pto_impedance = None
pto = wot.pto.PTO(3,kinematics, controller, pto_impedance, loss, name)  #      Here I set ndof = 3, and this block is run without error

#Constraints

def const_pos_pto_heave1(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps_p)[0]
    return 3 - np.abs(pos.flatten())

def const_pos_pto_heave2(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps_p)[1]
    return 3 - np.abs(pos.flatten())

def const_pos_pto_heave3(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps_p)[2]
    return 3 - np.abs(pos.flatten())

f_add = {'PTO': pto.force_on_wec}

constraints = [{'type': 'ineq','fun': const_pos_pto_heave1},
               {'type': 'ineq','fun': const_pos_pto_heave2},
               {'type': 'ineq','fun': const_pos_pto_heave3}]

#WEC

wec = wot.WEC.from_bem(
    bem_data,
    constraints=constraints,
    friction=None,
    f_add=f_add,
    inertia_matrix=[[40000,0,0],[0,40000,0],[0,0,40000]],
    hydrostatic_stiffness=[[100000,0,0],[0,10000,0],[0,0,10000]],
)

#Waves
amplitude = 0.75
wavefreq = 0.1
phase = 30
wavedir = 0
waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)

#Solve

obj_fun = pto.mechanical_average_power
nstate_opt = 2*nfreq

options = {'maxiter': 20}
scale_x_wec = 1e1
scale_x_opt = 1e-3
scale_obj = 1e-2

results = wec.solve(
    waves,
    obj_fun,
    nstate_opt,
    optim_options=options,
    scale_x_wec=scale_x_wec,
    scale_x_opt=scale_x_opt,
    scale_obj=scale_obj,
    )

opt_mechanical_average_power = results.fun
print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')

Error from email above:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[60], line 9
      6 scale_x_opt = 1e-3
      7 scale_obj = 1e-2
----> 9 results = wec.solve(
     10     waves,
     11     obj_fun,
     12     nstate_opt,
     13     optim_options=options,
     14     scale_x_wec=scale_x_wec,
     15     scale_x_opt=scale_x_opt,
     16     scale_obj=scale_obj,
     17     )
     19 opt_mechanical_average_power = results.fun
     20 print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')

File ~\anaconda3\envs\wot\Lib\site-packages\wecopttool\core.py:827, in WEC.solve(self, waves, obj_fun, nstate_opt, x_wec_0, x_opt_0, scale_x_wec, scale_x_opt, scale_obj, optim_options, use_grad, maximize, bounds_wec, bounds_opt, callback)
    824     problem['jac'] = grad(obj_fun_scaled)
    826 # minimize
--> 827 optim_res = minimize(**problem)
    829 msg = f'{optim_res.message}    (Exit mode {optim_res.status})'
    830 if optim_res.status == 0:

File ~\anaconda3\envs\wot\Lib\site-packages\scipy\optimize\_minimize.py:722, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    719     res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
    720                            bounds=bounds, **options)
    721 elif meth == 'slsqp':
--> 722     res = _minimize_slsqp(fun, x0, args, jac, bounds,
    723                           constraints, callback=callback, **options)
    724 elif meth == 'trust-constr':
    725     res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    726                                        bounds, constraints,
    727                                        callback=callback, **options)

File ~\anaconda3\envs\wot\Lib\site-packages\scipy\optimize\_slsqp_py.py:336, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
    322 exit_modes = {-1: "Gradient evaluation required (g & a)",
    323                0: "Optimization terminated successfully",
    324                1: "Function evaluation required (f & c)",
   (...)
    331                8: "Positive directional derivative for linesearch",
    332                9: "Iteration limit reached"}
    334 # Set the parameters that SLSQP will need
    335 # meq, mieq: number of equality and inequality constraints
--> 336 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    337           for c in cons['eq']]))
    338 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    339            for c in cons['ineq']]))
    340 # m = The total number of constraints

File ~\anaconda3\envs\wot\Lib\site-packages\scipy\optimize\_slsqp_py.py:336, in <listcomp>(.0)
    322 exit_modes = {-1: "Gradient evaluation required (g & a)",
    323                0: "Optimization terminated successfully",
    324                1: "Function evaluation required (f & c)",
   (...)
    331                8: "Positive directional derivative for linesearch",
    332                9: "Iteration limit reached"}
    334 # Set the parameters that SLSQP will need
    335 # meq, mieq: number of equality and inequality constraints
--> 336 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    337           for c in cons['eq']]))
    338 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    339            for c in cons['ineq']]))
    340 # m = The total number of constraints

File ~\anaconda3\envs\wot\Lib\site-packages\wecopttool\core.py:770, in WEC.solve.<locals>.scaled_resid_fun(x)
    768 x_s = x/scale
    769 x_wec, x_opt = self.decompose_state(x_s)
--> 770 return self._resid_fun(x_wec, x_opt, waves)

File ~\anaconda3\envs\wot\Lib\site-packages\wecopttool\core.py:611, in WEC._resid_fun(self, x_wec, x_opt, waves)
    609 # forces, -Σf
    610 for f in self.forces.values():
--> 611     ri = ri - f(self, x_wec, x_opt, waves)
    612 return self.dofmat_to_vec(ri)

File ~\anaconda3\envs\wot\Lib\site-packages\wecopttool\pto.py:341, in PTO.force_on_wec(self, wec, x_wec, x_opt, waves, nsubsteps)
    316 def force_on_wec(self,
    317     wec: TWEC,
    318     x_wec: ndarray,
   (...)
    321     nsubsteps: Optional[int] = 1,
    322 ) -> ndarray:
    323     """Calculate the PTO force on WEC.
    324
    325     Parameters
   (...)
    339         length.
    340     """
--> 341     force_td = self.force(wec, x_wec, x_opt, waves, nsubsteps)
    342     assert force_td.shape == (wec.nt*nsubsteps, self.ndof)
    343     force_td = np.expand_dims(np.transpose(force_td), axis=0)

File ~\anaconda3\envs\wot\Lib\site-packages\wecopttool\pto.py:124, in PTO.__init__.<locals>.force(wec, x_wec, x_opt, waves, nsubsteps)
    123 def force(wec, x_wec, x_opt, waves, nsubsteps=1):
--> 124     return controller(self, wec, x_wec, x_opt, waves, nsubsteps)

TypeError: 'list' object is not callable

[cmichelenstrofer]

So the biggest issue was with the shape of the BEM coefficients. They need to be ndof x ndof x nfreq. Everything becomes a 3x3 matrix. Capytaine has a neat function to create an array. The BEM is run for 3 DOFs (else there is not interaction and no need to model as an array). There were other small issues with shapes/indexing. This notebook works, simply replace in your geometry and mass.

For example, these are the inertia and hydrostatic stiffness matrices:

The constraints are respected:

Code:

import autograd.numpy as np
import wecopttool as wot

import capytaine as cpy
import matplotlib.pyplot as plt

# floating body

# single WEC
mesh = (wot.geom.WaveBot()).mesh(mesh_size_factor=0.5)
name = "WaveBot"
diameter = 1.76 # m

fb_base = cpy.FloatingBody.from_meshio(mesh, name=name)
fb_base.add_translation_dof(name="Heave")

# array
array = [
        [0, 0],
        [2 * diameter, 2 * diameter],
        [-2 * diameter, 2 * diameter],
    ]

fb = fb_base.assemble_arbitrary_array(np.array(array))
fb.name = "Array"
fb.mesh.name = "array"

ndof = fb.nb_dofs

#Bem_solver
f1 = 0.02
nfreq = 50
freq = wot.frequency(f1, nfreq, False)
bem_data = wot.run_bem(fb, freq)

# inertia_matrix
fb_base.center_of_mass = [0,0,0]
fb_base.rotation_center = fb_base.center_of_mass
mass = fb_base.compute_rigid_body_inertia().values[0][0]

bem_data.inertia_matrix[:] = np.diag([mass]*ndof)

bem_data

#PTO
name = ["PTO_Heave_0", "PTO_Heave_1", "PTO_Heave_2"]
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}

#Constraints
nsubsteps = 4

def const_pos_pto_heave1(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps)[:, 0]
    return 3 - np.abs(pos.flatten())

def const_pos_pto_heave2(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps)[:, 1]
    return 3 - np.abs(pos.flatten())

def const_pos_pto_heave3(wec, x_wec, x_opt, waves):
    pos = pto.position(wec, x_wec, x_opt, waves, nsubsteps)[:, 2]
    return 3 - np.abs(pos.flatten())

constraints = [{'type': 'ineq', 'fun': const_pos_pto_heave1},
               {'type': 'ineq', 'fun': const_pos_pto_heave2},
               {'type': 'ineq', 'fun': const_pos_pto_heave3}]

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

#Waves
amplitude = 0.75
wavefreq = 0.1
phase = 30
wavedir = 0
waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)

#Solve
obj_fun = pto.average_power
nstate_opt = 2*nfreq*ndof

options = {'maxiter': 100}
scale_x_wec = 1e1
scale_x_opt = 1e-3
scale_obj = 1e-2

results = wec.solve(
    waves,
    obj_fun,
    nstate_opt,
    optim_options=options,
    scale_x_wec=scale_x_wec,
    scale_x_opt=scale_x_opt,
    scale_obj=scale_obj,
    )

opt_mechanical_average_power = results[0].fun
print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')

# Analyze results
pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves.sel(realization=0), nsubsteps=nsubsteps)
wec_fdom, wec_tdom = wec.post_process(results[0], waves.sel(realization=0), nsubsteps=nsubsteps)

fig, ax = plt.subplots()
ax.plot(wec_tdom.time, np.ones(wec_tdom.time.shape)*3, 'k--')
ax.plot(wec_tdom.time, np.ones(wec_tdom.time.shape)*-3, 'k--')
wec_tdom.pos[0].plot(ax=ax, label="WEC_0")
wec_tdom.pos[1].plot(ax=ax, label="WEC_1")
wec_tdom.pos[2].plot(ax=ax, label="WEC_2")
ax.set_title("")
ax.set_ylabel("Heave [m]")
plt.legend()

Since the constraint is the same for all DOFs (+/- 3m) it could be simplified to a single constraint, but I left it as 3.