Kinematics (#63)

* make names of kinematics functions more explicit
* add a test for multiple DOF
* document the kinematics matrix

Co-authored-by: Ryan Coe <rcoe@sandia.gov>

CHANGES.md

@@ -4,12 +4,14 @@
 ## Version 1.0.3
 
 * minor updates to README
+* logging of decision vector and objective function
 * `f_add` should be passed as a `dict`, e.g., `{'my_name': my_func}`
   * optionally treat buoyancy/gravity explicitly via user-defined functions passed to `f_add`
   * time and freq domain results calculated for entries of `f_add` after `solve` completes
-* 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

@@ -109,6 +109,41 @@ The implicit assumption of this approach is that the body is neutrally buoyant (
 However, some WECs achieve equilibrium through a pretension applied via mooring and/or the PTO.
 In this case, the device can still be modeled with the linear hydrostatic stiffness, but if you wish to solve for the pretension force in your simulations, you may explicitly include the buoyancy, gravity, and pretension forces via the :code:`f_add` argument to :py:class:`wecopttool.core.WEC`.
 
+PTO Kinematics
+^^^^^^^^^^^^^^
+The :py:mod:`wecopttool.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*.
+Here, the kinematics matrix, :math:`K`, is defined as the linear transformation from the WEC position (e.g., heave) in the global frame, :math:`x`, to the PTO position in the PTO frame (e.g., tether length/generator rotation), :math:`p`:
+
+.. math::
+    p = K x
+    :label: kinematics
+
+The relationship :math:`p(x)` is typically referred to as the *forward kinematics*.
+The matrix :math:`K` has a size equal to the number of DOFs of the PTOs times the number of DOFs of the WEC.
+Note, however that the real kinematics might not be linear.
+Equation :eq:`kinematics` represents a linearization of :math:`p(x)` about the mean :math:`x=0` position, with the matrix :math:`K` being the Jacobian of :math:`p(x)` at :math:`x=0`.
+
+The transpose of :math:`K` is used to transform the PTO forces in PTO frame, :math:`f_p`, to the PTO forces on the WEC, :math:`f_w`:
+
+.. math::
+    f_w = K^T f_p
+    :label: k
+
+This relationship can be derived from conservation of energy in both frames, and using the definition in Equation :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
+
+..
+    (commented out): This represents a linearization of the function :math:`f_w(f_p)` about :math:`f_p=0` with :math:`K^T` being the Jacobian of :math:`f_w(f_p)` at :math:`f_p=0`.
+                     The assumption here is that :math:`f_p(p(x=0))=f_p(0)=0`.
+
+
 .. _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()
 
@@ -322,32 +322,159 @@ def test_examples_device_wavebot_plot_cross_section():
     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': 250, 'ftol': 1e-8}
+    scale_x_wec = 100.0
+    scale_x_opt = 0.1
+    scale_obj = 1.0
+    _, _, _, _, avg_pow_sh_sh, _ = wec.solve(
+        regular_wave, obj_fun, nstate_opt, optim_options=options,
+        scale_x_wec=scale_x_wec, scale_x_opt=scale_x_opt, scale_obj=scale_obj)
+
+    # 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,
+        scale_x_wec=scale_x_wec, scale_x_opt=scale_x_opt, scale_obj=scale_obj)
+
+    # 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,
+        scale_x_wec=scale_x_wec, scale_x_opt=scale_x_opt, scale_obj=scale_obj)
+
+    # 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,
+        scale_x_wec=scale_x_wec, scale_x_opt=scale_x_opt, scale_obj=scale_obj)
+
+    # 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
+    scale_x_wec = 10000.0
+    scale_x_opt = 1.0
+    scale_obj = 1.0
+    _, _, _, _, avg_pow_h_h, _ = wec.solve(
+        regular_wave, obj_fun, nstate_opt, optim_options=options,
+        scale_x_wec=scale_x_wec, scale_x_opt=scale_x_opt, scale_obj=scale_obj)
+
+    # checks
+    tol = 1e0
+    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
+
+
 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"""
-    
+
     delta = np.random.randn() # excess buoyancy factor
-    
+
     # remove constraints
     wec.constraints = []
-    
+
     def f_b(wec, x_wec, x_opt):
         V = wec.fb.keep_immersed_part(inplace=False).mesh.volume
         rho = wec.rho
         g = 9.81
         return (1+delta) * rho * V * g * np.ones([wec.ncomponents, wec.ndof])
-    
+
     def f_g(wec, x_wec, x_opt):
         g = 9.81
         m = wec.mass.item()
         return -1 * m * g * np.ones([wec.ncomponents, wec.ndof])
-    
+
     wec.f_add = {**wec.f_add,
                  'Fb':f_b,
                  'Fg':f_g,
                  }
-    
-    tdom, fdom, x_wec, x_opt, avg_pow, _ = wec.solve(regular_wave, 
+
+    tdom, fdom, x_wec, x_opt, avg_pow, _ = wec.solve(regular_wave,
                                                 obj_fun = pto.average_power,
                                                 nstate_opt = pto.nstate,
                                                 scale_x_wec = 1.0,
@@ -359,4 +486,3 @@ def f_g(wec, x_wec, x_opt):
     expected = (wec.rho * wec.fb.mesh.volume * wec.g * delta) \
         / wec.hydrostatic_stiffness.item()
     assert pytest.approx (expected, 1e-1) == mean_pos
-

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").