Open En-Lo opened 9 months ago
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()
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:
import autograd.numpy as np
instead of just numpyOther issues:
Recommendations:
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
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 #
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_tdvelpto_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
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, ...]
.
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:
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:
To Reproduce
System: