Initial guess, automatic scaling, multi-DOF power limit

tests/test_wecopttool.py

@@ -137,6 +137,36 @@ def test_solve(wec, regular_wave, pto):
     _, _ = pto.post_process(wec, x_wec, x_opt)
 
 
+def test_solve_initial_guess(wec, regular_wave, pto):
+
+    x_wec_0 = wec.initial_x_wec_guess(regular_wave)
+
+    nits = []
+    x_wecs = []
+    for x_wec_0i in [None, x_wec_0]:
+        *_, x_wec, _, _, res = wec.solve(regular_wave,
+                            obj_fun=pto.average_power,
+                            nstate_opt=pto.nstate,
+                            scale_x_wec=1.0,
+                            scale_x_opt=0.01,
+                            scale_obj=1e-1,
+                            x_wec_0=x_wec_0i,
+                            )
+        nits.append(res['nit'])
+        x_wecs.append(x_wec)
+
+    assert nits[0] > nits[1]
+    assert pytest.approx(x_wecs[1],1e0) == x_wec_0
+
+
+def test_complex_to_real_amplitudes(wec, regular_wave):
+    x_wec = wec.initial_x_wec_guess(regular_wave)
+    fd_wec = wot.real_to_complex_amplitudes(np.atleast_2d(x_wec))
+    x_wec_1 = wot.complex_to_real_amplitudes(fd_wec)
+
+    assert np.all(x_wec == x_wec_1)
+
+
 def test_solve_constraints(wec, regular_wave, pto):
     """Checks that two constraints on PTO force can be enforced
     """
@@ -260,9 +290,15 @@ def test_wavebot_ps_theoretical_limit(wec,regular_wave,pto):
     nstate_opt = pto.nstate
     _, fdom, _, _, average_power, _ = wec.solve(regular_wave, obj_fun, nstate_opt,
         optim_options={'maxiter': 1000, 'ftol': 1e-8}, scale_x_opt=1e3)
-    plim = power_limit(fdom['excitation_force'][1:, 0], wec.hydro.Zi[:, 0, 0])
-
+
+    plim = power_limit(fdom['excitation_force'], wec.hydro.Zi)
     assert pytest.approx(average_power, 1e-5) == plim
+
+    opt_vel = wec.optimal_velocity(regular_wave)
+    assert pytest.approx(fdom.vel.data[1:],1e0) == opt_vel
+
+    opt_pos = wec.optimal_position(regular_wave)
+    assert pytest.approx(fdom.pos.data[1:],1e0) == opt_pos
 
 
 def test_wavebot_p_cc(wec,resonant_wave):
@@ -289,8 +325,8 @@ def test_wavebot_p_cc(wec,resonant_wave):
         bounds=bounds)
 
     # P controller power matches theoretical limit at resonance
-    plim = power_limit(fdom['excitation_force'][1:, 0],
-                       wec.hydro.Zi[:, 0, 0]).item()
+    plim = power_limit(fdom['excitation_force'],
+                       wec.hydro.Zi).item()
 
     assert pytest.approx(average_power, 0.03) == plim
 
@@ -331,8 +367,8 @@ def test_wavebot_pi_cc(wec,regular_wave):
                                                bounds=bounds)
 
     # PI controller power matches theoretical limit for an single freq
-    plim = power_limit(fdom['excitation_force'][1:, 0],
-                       wec.hydro.Zi[:, 0, 0]).item()
+    plim = power_limit(fdom['excitation_force'],
+                       wec.hydro.Zi).item()
 
     assert pytest.approx(avg_power, 0.03) == plim
 

wecopttool/core.py

@@ -8,6 +8,7 @@
 __all__ = ['WEC', 'freq_array', 'real_to_complex_amplitudes', 'fd_to_td',
            'td_to_fd', 'scale_dofs', 'complex_xarray_from_netcdf',
            'complex_xarray_to_netcdf', 'wave_excitation', 'run_bem',
+           'optimal_velocity', 'optimal_position', 'complex_to_real_amplitudes',
            'power_limit', 'natural_frequency', 'plot_impedance',
            'post_process_continuous_time']
 
@@ -579,6 +580,17 @@ def plot_impedance(self, style: str = 'Bode', option: str = 'symmetric',
             show=show)
         return fig, axs
 
+    def optimal_velocity(self, waves: xr.DataSet):
+        fd_wec, _ = wave_excitation(self.hydro, waves)
+        return optimal_velocity(excitation=fd_wec['excitation_force'],
+                                impedance=self.hydro['Zi'])
+
+    def optimal_position(self, waves: xr.DataSet):
+        fd_wec, _ = wave_excitation(self.hydro, waves)
+        return optimal_position(excitation=fd_wec['excitation_force'],
+                                impedance=self.hydro['Zi'],
+                                omega=self.hydro['omega'])
+
     def natural_frequency(self):
         """Return natural frequency or frequencies.
 
@@ -608,6 +620,25 @@ def _get_state_scale(self,
 
         return np.concatenate([scale_x_wec, scale_x_opt])
 
+    def initial_x_wec_guess(self, waves: xr.Dataset) -> np.ndaray:
+        """Initial guess for `x_wec` based on optimal hydrodynamic solution to 
+        be passed to `wec.solve`.
+
+        Parameters
+        ----------
+        waves : xr.Dataset
+            Wave DataSet
+
+        Returns
+        -------
+        x_wec_0
+            Initial guess for `x_wec`
+        """
+        pos_opt = self.optimal_position(waves)
+        pos_opt_zero_mean = np.concatenate([np.zeros((1, self.ndof)), pos_opt])
+        x_wec_0 = complex_to_real_amplitudes(pos_opt_zero_mean)
+        return x_wec_0.squeeze()
+
     def solve(self,
               waves: xr.Dataset,
               obj_fun: Callable[[WEC, np.ndarray, np.ndarray], float],
@@ -904,6 +935,21 @@ def real_to_complex_amplitudes(fd: np.ndarray, first_row_is_mean: bool = True
     return np.concatenate((mean, fd[0::2, :] - 1j*fd[1::2, :]), axis=0)
 
 
+def complex_to_real_amplitudes(fd: np.ndarray) -> np.ndarray:
+    """Convert from one complex amplitude to two real amplitudes per 
+    frequency."""
+
+    m = fd.shape[0]
+    n = fd.shape[1]
+    out = np.zeros((1+2*(m-1),n))
+
+    out[0,:] = fd[0,:]
+    out[1::2,:] = fd[1:].real
+    out[2::2,:] = -1*fd[1:].imag
+
+    return out
+
+
 def fd_to_td(fd: np.ndarray, n: Optional[int] = None) -> np.ndarray:
     return np.fft.irfft(fd/2, n=n, axis=0, norm='forward')
 
@@ -934,7 +980,7 @@ def scale_dofs(scale_list: list[float], ncomponents: int) -> np.ndarray:
 
 
 def complex_xarray_from_netcdf(fpath: str | Path) -> xr.Dataset:
-    """Read a NetCDF file with complex entries as an xarray dataSet.
+    """Read a NetCDF file with complex entries as an xarray DataSet.
     """
     with xr.open_dataset(fpath) as ds:
         ds.load()
@@ -1088,14 +1134,57 @@ def run_bem(fb: cpy.FloatingBody, freq: Iterable[float] = [np.infty],
     return solver.fill_dataset(test_matrix, [wec_im], **write_info)
 
 
+def optimal_velocity(excitation: npt.ArrayLike, impedance: npt.ArrayLike
+                     ) -> np.ndarray:
+    """Find optimal velocity.
+
+    Parameters
+    ----------
+    excitation: np.ndarray
+        Complex excitation spectrum. Shape: ``nfreq`` x ``ndof``
+    impedance: np.ndarray
+        Complex impedance matrix. Shape: ``nfreq`` x ``ndof`` x ``ndof``
+
+    Returns
+    -------
+    opt_vel
+        Optimal velocity for power absorption.
+    """
+    opt_vel = np.concatenate([np.linalg.lstsq(2*impedance[w_ind, :, :].real,
+                                              excitation[w_ind+1, :])[0] 
+                              for w_ind in range(impedance.shape[0])])
+    return np.atleast_2d(opt_vel).transpose()
+
+
+def optimal_position(excitation: npt.ArrayLike, impedance: npt.ArrayLike,
+                     omega: npt.ArrayLike) -> np.ndarray:
+    """Find optimal position.
+
+    Parameters
+    ----------
+    excitation: np.ndarray
+        Complex excitation spectrum. Shape: ``nfreq`` x ``ndof``
+    impedance: np.ndarray
+        Complex impedance matrix. Shape: ``nfreq`` x ``ndof`` x ``ndof``
+
+    Returns
+    -------
+    optimal_position
+        Optimal position for power absorption.
+    """
+    opt_vel = optimal_velocity(excitation, impedance)
+    opt_pos = opt_vel / (1j * np.atleast_2d(omega.data).transpose())
+    return opt_pos
+
+
 def power_limit(excitation: npt.ArrayLike, impedance: npt.ArrayLike
                 ) -> np.ndarray:
     """Find upper limit for power.
 
     Parameters
     ----------
-    exctiation: np.ndarray
-        Complex exctitation spectrum. Shape: ``nfreq`` x ``ndof``
+    excitation: np.ndarray
+        Complex excitation spectrum. Shape: ``nfreq`` x ``ndof``
     impedance: np.ndarray
         Complex impedance matrix. Shape: ``nfreq`` x ``ndof`` x ``ndof``
 
@@ -1106,6 +1195,12 @@ def power_limit(excitation: npt.ArrayLike, impedance: npt.ArrayLike
     """
 
     power_limit = -1*np.sum(np.abs(excitation)**2 / (8*np.real(impedance)))
+
+    pls = np.concatenate([np.linalg.lstsq(8*impedance[w_ind, :, :].real,
+                                          np.abs(excitation[w_ind+1, :])**2)[0] 
+                         for w_ind in range(impedance.shape[0])])
+
+    power_limit = -1 * np.sum(pls)
 
     return power_limit