Kinematics

CHANGES.md

@@ -4,9 +4,10 @@
 ## Version 1.0.3
 
 * minor updates to README
-* logging of decision vector and objective function 
+* logging of decision vector and objective function
   * 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.
 
 **Bug fixes**
 

docs/source/theory.rst

@@ -102,6 +102,33 @@ An important practical factor when using this functionality is to make sure that
 It may be required to enforce constraints at more points than the dynamics (as defined by the frequency array).
 In WecOptTool's example PTO module, this is controlled by the :code:`nsubsteps` argument (see, e.g., :meth:`pto.PseudoSpectralPTO.force`).
 
+PTO Kinematics
+^^^^^^^^^^^^^^
+The PTO module includes several examples of PTOs that can be used for both additional PTO forces on the WEC dynamics and for objective functions (e.g., PTO average power).
+Creating one of these pre-defined PTOs requires specifying the *kinematics matrix*.
+The kinematics matrix :math:`K` transforms the WEC position :math:`x` (e.g. heave) in the global frame to the PTO position :math:`p` in the PTO frame (e.g. tether length/generator rotation).
+It has size equal to the number of DOFs of the PTOs times the number of DOFs of the WEC.
+It is used to obtain the PTO motion as:
+
+.. math::
+    p = K x
+    :label: kinematics
+
+The kinematics matrix is the transpose of the Jacobian matrix used to transform the PTO force :math:`f_p` in PTO frame to the PTO forces on the WEC :math:`f_w` in the global frame, as
+
+.. math::
+    f_w = K^T f_p
+    :label: k
+
+This relation can be derived from conservation of energy in both frames, and using :eq:`kinematics`:
+
+.. math::
+    f_w^T x = f_p^T p \\
+    f_w^T x = f_p^T K x \\
+    f_w^T = f_p^T K \\
+    f_w = K^T f_p \\
+    :label: conservation_energy
+
 .. _WEC-Sim: https://wec-sim.github.io/WEC-Sim/master/index.html
 .. _Autograd: https://github.com/HIPS/autograd
 .. _Autograd documentation: https://github.com/HIPS/autograd/blob/master/docs/tutorial.md#supported-and-unsupported-parts-of-numpyscipy

tests/test_wecopttool.py

@@ -139,7 +139,7 @@ def test_solve(wec, regular_wave, pto):
 def test_solve_constraints(wec, regular_wave, pto):
     """Checks that two constraints on PTO force can be enforced
     """
-    
+
     f_max = 1.85e3
     f_min = 1.8e3
     nsubsteps = 4
@@ -151,17 +151,17 @@ def const_f_pto_max(wec, x_wec, x_opt):
     def const_f_pto_min(wec, x_wec, x_opt):
         f = pto.force_on_wec(wec, x_wec, x_opt, nsubsteps)
         return f.flatten() + f_min
-        
+
     ineq_cons_max = {'type': 'ineq',
                 'fun': const_f_pto_max,
                 }
 
     ineq_cons_min = {'type': 'ineq',
                 'fun': const_f_pto_min,
                 }
-    
+
     wec.constraints = [ineq_cons_max, ineq_cons_min]
-    
+
     options = {}
     obj_fun = pto.average_power
     nstate_opt = pto.nstate
@@ -171,9 +171,9 @@ def const_f_pto_min(wec, x_wec, x_opt):
         scale_x_opt = 0.01,
         scale_obj = 1e-1,
         optim_options=options)
-    
+
     pto_tdom, _ = pto.post_process(wec, x_wec, x_opt)
-    
+
     assert pytest.approx(-1*f_min, 1e-5) == pto_tdom['force'].min().values.item()
     assert pytest.approx(f_max, 1e-5) == pto_tdom['force'].max().values.item()
 
@@ -320,3 +320,119 @@ def test_examples_device_wavebot_mesh():
 def test_examples_device_wavebot_plot_cross_section():
     wb = WaveBot()
     wb.plot_cross_section()
+
+
+def test_multiple_dof(regular_wave):
+    """When defined as uncoupled, surge and heave computed seperately should 
+    give same solution as computed together"""
+
+    # 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]])
+
+    # PTO
+    kinematics = np.eye(fb.nb_dofs)
+    names = ["SURGE", "HEAVE"]
+    pto = wot.pto.PseudoSpectralPTO(nfreq, kinematics, names=names)
+
+    # WEC
+    f_add = pto.force_on_wec
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq, f_add=f_add)
+
+    # 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()
+
+    # solve
+    obj_fun = pto.average_power
+    nstate_opt = pto.nstate
+    options = {'maxiter': 1000, 'ftol': 1e-8}
+    _, _, _, _, avg_pow_sh_sh, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
+                                             optim_options=options)
+
+    # only surge PTO
+    kinematics = np.array([[1.0, 0.0]])
+    pto = wot.pto.PseudoSpectralPTO(nfreq, kinematics, names=["SURGE"])
+    f_add = pto.force_on_wec
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq, f_add=f_add)
+    wec.run_bem()
+    obj_fun = pto.average_power
+    nstate_opt = pto.nstate
+    _, _, _, _, avg_pow_sh_s, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
+                                            optim_options=options)
+
+    # only heave PTO
+    kinematics = np.array([[0.0, 1.0]])
+    pto = wot.pto.PseudoSpectralPTO(nfreq, kinematics, names=["HEAVE"])
+    f_add = pto.force_on_wec
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq, f_add=f_add)
+    wec.run_bem()
+    obj_fun = pto.average_power
+    nstate_opt = pto.nstate
+    _, _, _, _, avg_pow_sh_h, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
+                                            optim_options=options)
+
+    # WEC only surge
+    mass = np.array([[mass_11]])
+    stiffness = np.array([[stiffness_11]])
+    fb = cpy.FloatingBody.from_file(meshfile, name="WaveBot")
+    fb.add_translation_dof(name="SURGE")
+    kinematics = np.array([[1.0]])
+    pto = wot.pto.PseudoSpectralPTO(nfreq, kinematics, names=["SURGE"])
+    f_add = pto.force_on_wec
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq, f_add=f_add)
+    wec.run_bem()
+    obj_fun = pto.average_power
+    nstate_opt = pto.nstate
+    _, _, _, _, avg_pow_s_s, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
+                                           optim_options=options)
+
+    # WEC only heave
+    mass = np.array([[mass_33]])
+    stiffness = np.array([[stiffness_33]])
+    fb = cpy.FloatingBody.from_file(meshfile, name="WaveBot")
+    fb.add_translation_dof(name="HEAVE")
+    pto = wot.pto.PseudoSpectralPTO(nfreq, kinematics, names=["HEAVE"])
+    f_add = pto.force_on_wec
+    wec = wot.WEC(fb, mass, stiffness, f0, nfreq, f_add=f_add)
+    wec.run_bem()
+    obj_fun = pto.average_power
+    nstate_opt = pto.nstate
+    _, _, _, _, avg_pow_h_h, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
+                                           optim_options=options)
+
+    # checks
+    tol = 1e-2
+    assert pytest.approx(avg_pow_sh_sh, tol) == avg_pow_sh_s + avg_pow_sh_h
+    assert pytest.approx(avg_pow_sh_s, tol) == avg_pow_s_s
+    assert pytest.approx(avg_pow_sh_h, tol) == avg_pow_h_h

wecopttool/pto.py

@@ -60,26 +60,26 @@ def nstate(self):
     def _kinematics_t(self):
         return np.transpose(self.kinematics)
 
-    def _pto_to_wec_dofs(self, pto):
-        return np.dot(pto, self.kinematics)
+    def _pto_force_to_wec_force(self, f_pto):
+        return np.dot(f_pto, self.kinematics)
 
-    def _wec_to_pto_dofs(self, pto):
-        return np.dot(pto, self._kinematics_t)
+    def _wec_pos_to_pto_pos(self, pos_wec):
+        return np.dot(pos_wec, self._kinematics_t)
 
     def position(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
                  nsubsteps: int = 1) -> np.ndarray:
         """Calculate the PTO position time-series."""
         wec_pos = wec.vec_to_dofmat(x_wec)
-        pos = self._wec_to_pto_dofs(wec_pos)
+        pos = self._wec_pos_to_pto_pos(wec_pos)
         time_mat = wec.make_time_mat(nsubsteps)
         return np.dot(time_mat, pos)
 
     def velocity(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
                  nsubsteps: int = 1) -> np.ndarray:
         """Calculate the PTO velocity time-series."""
         wec_pos = wec.vec_to_dofmat(x_wec)
-        wec_vel = np.dot(wec.derivative_mat, wec_pos)
-        vel = self._wec_to_pto_dofs(wec_vel)
+        pos = self._wec_pos_to_pto_pos(wec_pos)
+        vel = np.dot(wec.derivative_mat, pos)
         time_mat = wec.make_time_mat(nsubsteps)
         return np.dot(time_mat, vel)
 
@@ -88,9 +88,9 @@ def acceleration(self, wec: WEC, x_wec: npt.ArrayLike,
                      ) -> np.ndarray:
         """Calculate the PTO acceleration time-series."""
         wec_pos = wec.vec_to_dofmat(x_wec)
-        wec_vel = np.dot(wec.derivative_mat, wec_pos)
-        wec_acc = np.dot(wec.derivative_mat, wec_vel)
-        acc = self._wec_to_pto_dofs(wec_acc)
+        pos = self._wec_pos_to_pto_pos(wec_pos)
+        vel = np.dot(wec.derivative_mat, pos)
+        acc = np.dot(wec.derivative_mat, vel)
         time_mat = wec.make_time_mat(nsubsteps)
         return np.dot(time_mat, acc)
 
@@ -111,7 +111,7 @@ def force_on_wec(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
         See ``force``.
         """
         fpto_td = self.force(wec, x_wec, x_opt, nsubsteps)
-        return self._pto_to_wec_dofs(fpto_td)
+        return self._pto_force_to_wec_force(fpto_td)
 
     def energy(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
                nsubsteps: int = 1) -> float:
@@ -302,15 +302,16 @@ def energy(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
                nsubsteps: int = 1) -> float:
         if nsubsteps == 1:
             wec_pos = wec.vec_to_dofmat(x_wec)
-            wec_vel = np.dot(wec.derivative_mat, wec_pos)
-            vel = self._wec_to_pto_dofs(wec_vel)
+            pos = self._wec_pos_to_pto_pos(wec_pos)
+            vel = np.dot(wec.derivative_mat, pos)
             vel_vec = wec.dofmat_to_vec(vel[1:, :])
             energy_produced = 1/(2*wec.f0) * np.dot(vel_vec, x_opt)
         else:
             energy_produced = super().energy(wec, x_wec, x_opt, nsubsteps)
         return energy_produced
 
 
+
 class ProportionalPTO(_PTO):
     """Proportional (P) PTO controller (a.k.a., "proportional damping").