BUG fix: structured controllers (#71)

structured controller (PI and P) PTOs were not working properly for multi-DOF WECs

CHANGES.md

@@ -1,7 +1,7 @@
 
 # Changelog
 
-## Version 1.0.3
+## Version 1.1.0
 
 * minor updates to README
 * logging of decision vector and objective function
@@ -12,11 +12,13 @@
   * controlled entirely via logging package config
   * move to `info` logging level (`debug` gives too much from other packages, e.g., matplotlib)
 * added tests for multiple WEC/PTO degrees of freedom.
+* allow user to pass `bounds` via `solve` to `scipy.optimize.minimize`
 
 **Bug fixes**
 
 * geom.WaveBot.plot_cross_section now plots correct location of waterline
 * All constraints are now being enforced when multiple constraints are passed
+* Fix shape/linear algebra bugs in fixed structure PTOs with multiple PTO DOFs
 
 
 ## Version 1.0.2

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = wecopttool
-version = 1.0.2
+version = 1.1.0
 author = Sandia National Laboratories
 description = WEC Design Optimization Toolbox
 long_description = file: README.md

tests/test_wecopttool.py

@@ -8,6 +8,7 @@
 from autograd.builtins import isinstance, tuple, list, dict
 import capytaine as cpy
 import meshio
+from scipy.optimize import Bounds
 
 import wecopttool as wot
 from wecopttool.geom import WaveBot
@@ -279,8 +280,16 @@ def test_wavebot_p_cc(wec,resonant_wave):
     nstate_opt = pto.nstate
     wec.f_add = {'PTO': pto.force_on_wec}
 
+    # set bounds such that damping must be negative
+    lb = np.concatenate([-1 * np.inf * np.ones(wec.nstate_wec), 
+                         -1 * np.inf * np.ones(nstate_opt)])
+    ub = np.concatenate([1 * np.inf * np.ones(wec.nstate_wec), 
+                         0 * np.ones(nstate_opt)])
+    bounds = Bounds(lb, ub)
+
     _, fdom, _, xopt, average_power, _ = wec.solve(resonant_wave, obj_fun, nstate_opt,
-        optim_options={'maxiter': 1000, 'ftol': 1e-8}, scale_x_opt=1e3)
+        optim_options={'maxiter': 1000, 'ftol': 1e-8}, scale_x_opt=1e3,
+        bounds=bounds)
     plim = power_limit(fdom['excitation_force'][1:, 0],
                        wec.hydro.Zi[:, 0, 0]).item()
 
@@ -449,6 +458,87 @@ def test_multiple_dof(regular_wave):
     assert pytest.approx(avg_pow_sh_h, tol) == avg_pow_h_h
 
 
+@pytest.fixture
+def surge_heave_wavebot():
+    # frequencies
+    f0 = 0.05
+    nfreq = 18
+
+    #  mesh
+    meshfile = os.path.join(os.path.dirname(__file__), 'data', 'wavebot.stl')
+
+    # capytaine floating body
+    fb = cpy.FloatingBody.from_file(meshfile, name="WaveBot")
+    fb.add_translation_dof(name="SURGE")
+    fb.add_translation_dof(name="HEAVE")
+
+    # hydrostatic
+    hs_data = wot.hydrostatics.hydrostatics(fb)
+    mass_11 = wot.hydrostatics.mass_matrix_constant_density(hs_data)[0, 0]
+    mass_13 = wot.hydrostatics.mass_matrix_constant_density(hs_data)[0, 2] # will be 0
+    mass_31 = wot.hydrostatics.mass_matrix_constant_density(hs_data)[2, 0] # will be 0
+    mass_33 = wot.hydrostatics.mass_matrix_constant_density(hs_data)[2, 2]
+    stiffness_11 = wot.hydrostatics.stiffness_matrix(hs_data)[0, 0]
+    stiffness_13 = wot.hydrostatics.stiffness_matrix(hs_data)[0, 2] # will be 0
+    stiffness_31 = wot.hydrostatics.stiffness_matrix(hs_data)[2, 0] # will be 0
+    stiffness_33 = wot.hydrostatics.stiffness_matrix(hs_data)[2, 2]
+    mass = np.array([[mass_11, mass_13],
+                     [mass_31, mass_33]])
+    stiffness = np.array([[stiffness_11, stiffness_13],
+                          [stiffness_31, stiffness_33]])
+
+    # WEC
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq)
+
+    # BEM
+    wec.run_bem()
+
+    # set diagonal (coupling) components to zero
+    wec.hydro.added_mass[:, 0, 1] = 0.0
+    wec.hydro.added_mass[:, 1, 0] = 0.0
+    wec.hydro.radiation_damping[:, 0, 1] = 0.0
+    wec.hydro.radiation_damping[:, 1, 0] = 0.0
+    wec._del_impedance()
+    wec.bem_calc_impedance()
+
+    return wec
+
+
+def test_multiple_dof_fixed_structure_P(regular_wave, surge_heave_wavebot):
+
+    kinematics = np.eye(surge_heave_wavebot.ndof)
+    names = ["SURGE", "HEAVE"]
+    pto = wot.pto.ProportionalPTO(kinematics, names=names)
+
+    x_wec = np.random.randn(surge_heave_wavebot.nstate_wec)
+    x_opt = np.random.randn(pto.nstate)
+
+    pto_force = pto.force_on_wec(surge_heave_wavebot, x_wec, x_opt)
+    pto_vel = pto.velocity(surge_heave_wavebot, x_wec, x_opt)
+
+    assert np.all(x_opt[0] * pto_vel[:,0] == pto_force[:,0])
+    assert np.all(x_opt[1] * pto_vel[:,1] == pto_force[:,1])
+
+
+def test_multiple_dof_fixed_structure_PI(regular_wave, surge_heave_wavebot):
+
+    kinematics = np.eye(surge_heave_wavebot.ndof)
+    names = ["SURGE", "HEAVE"]
+    pto = wot.pto.ProportionalIntegralPTO(kinematics, names=names)
+
+    x_wec = np.random.randn(surge_heave_wavebot.nstate_wec)
+    x_opt = np.random.randn(pto.nstate)
+
+    pto_force = pto.force_on_wec(surge_heave_wavebot, x_wec, x_opt)
+    pto_vel = pto.velocity(surge_heave_wavebot, x_wec, x_opt)
+    pto_pos = pto.position(surge_heave_wavebot, x_wec, x_opt)
+
+    assert np.all(x_opt[0] * pto_vel[:,0] + x_opt[2] * pto_pos[:,0] 
+                  == pto_force[:,0])
+    assert np.all(x_opt[1] * pto_vel[:,1] + x_opt[3] * pto_pos[:,1] 
+                  == pto_force[:,1])
+
+
 def test_buoyancy_excess(wec, pto, regular_wave):
     """Give too much buoyancy and check that equilibrium point found matches
     that given by the hydrostatic stiffness"""

wecopttool/core.py

@@ -22,7 +22,7 @@
 from autograd import grad, jacobian
 import xarray as xr
 import capytaine as cpy
-from scipy.optimize import minimize, OptimizeResult
+from scipy.optimize import minimize, OptimizeResult, Bounds
 from scipy.linalg import block_diag
 import matplotlib.pyplot as plt
 import matplotlib as mpl
@@ -620,6 +620,7 @@ def solve(self,
               optim_options: dict[str, Any] = {},
               use_grad: bool = True,
               maximize: bool = False,
+              bounds: Optional[Bounds | list] = None,
               ) -> tuple[xr.Dataset, xr.Dataset, np.ndarray, np.ndarray, float,
                          OptimizeResult]:
         """Solve the WEC co-design problem.
@@ -665,6 +666,8 @@ def solve(self,
         maximize: bool
             Whether to maximize the objective function. The default is
             ``False`` to minimize the objective function.
+        bounds: sequence | Bounds
+            See scipy.optimize.minimize
 
         Returns
         -------
@@ -741,6 +744,7 @@ def resid_fun(x):
                    'method': 'SLSQP',
                    'constraints': constraints,
                    'options': optim_options,
+                   'bounds': bounds,
                    }
 
         def callback(x):

wecopttool/pto.py

@@ -332,7 +332,7 @@ def nstate_per_dof(self):
     def force(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
               nsubsteps: int = 1) -> np.ndarray:
         vel_td = self.velocity(wec, x_wec, x_opt, nsubsteps)
-        force_td = np.reshape(x_opt, [-1, 1]) * vel_td
+        force_td = np.reshape(x_opt, [1,-1]) * vel_td
         return force_td
 
 
@@ -355,7 +355,8 @@ def force(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
               nsubsteps: int = 1) -> np.ndarray:
         vel_td = self.velocity(wec, x_wec, x_opt, nsubsteps)
         pos_td = self.position(wec, x_wec, x_opt, nsubsteps)
-        u = np.reshape(x_opt, [-1, 1])
+        u = np.reshape(x_opt, [1,-1])
         B = np.hstack([vel_td, pos_td])
-        force_td = np.dot(B,u)
+        tmp1 = u * B
+        force_td = tmp1[:,0:self.ndof] + tmp1[:,self.ndof:]
         return force_td