diff --git a/docs/source/theory.rst b/docs/source/theory.rst
index c5b9f59cf..c8169c576 100644
--- a/docs/source/theory.rst
+++ b/docs/source/theory.rst
@@ -146,6 +146,16 @@ This relationship can be derived from conservation of energy in both frames, and
f_w = K^T f_p \\
:label: conservation_energy
+Phase Realizations
+^^^^^^^^^^^^^^^^^^
+Irregular waves are defined in WecOptTool as a spectrum of complex frequency-domain wave elevations. The phase of each of the elevation elements impacts the time-domain result. Thus, the standard calculation of the objective function (average power) may change across a range of phase realizations. The amount of variance in power depends on multiple factors such as the optimization problem and the frequency array. When creating an irregular wave using :py:meth:`wecopttool.waves.long_crested_wave` or :py:meth:`wecopttool.waves.irregular_wave`, :code:`nrealizations` can be used to select the number of phase realizations to be used. By default, random realizations will be used to create the selected number of wave elevation spectra. The :py:meth:`wecopttool.WEC.solve` function will automatically iterate through and solve the optimization problem for each realization, and the overall result can be found by averaging the value of the objective function across all realizations. A general recommendation is to use sufficient random phase realizations such that the total simulation time sums to 20 minutes.
+
+.. math::
+ t_{total} = \frac{nrealizations}{f1}
+ :label: total_time
+
+The selection of the number of realizations is further detailed in :doc:`_examples/tutorial_4_Pioneer`.
+
Troubleshooting
---------------
If your simulation is not behaving as expected, consider some of the general troubleshooting steps below:
diff --git a/examples/data/tutorial_4_freq_range.nc b/examples/data/tutorial_4_freq_range.nc
new file mode 100644
index 000000000..37a976e71
Binary files /dev/null and b/examples/data/tutorial_4_freq_range.nc differ
diff --git a/examples/data/tutorial_4_nfreqs.nc b/examples/data/tutorial_4_nfreqs.nc
new file mode 100644
index 000000000..c2c9c5d13
Binary files /dev/null and b/examples/data/tutorial_4_nfreqs.nc differ
diff --git a/examples/data/tutorial_4_nrealizations.nc b/examples/data/tutorial_4_nrealizations.nc
new file mode 100644
index 000000000..b0029b696
Binary files /dev/null and b/examples/data/tutorial_4_nrealizations.nc differ
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index 030077bab..7be4ff5aa 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -115,12 +115,43 @@
"_ = wb.plot_cross_section(show=True) # specific to WaveBot"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Waves\n",
+ "In this case we will use a regular wave with a frequency of 0.3 Hz and an amplitude of 0.0625 m. \n",
+ "\n",
+ "In regular waves, the WEC system response will occur at the wave frequency and corresponding odd harmonics (3rd, 5th, etc.) of the wave frequency (due to nonlinearities such as PTO force constraints). Thus, we can set the fundamental frequency, $f_1$, (also equal to the array spacing) equal to the wave frequency. We analyzed the system response in the frequency domain and determined that 10 frequencies was sufficient to capture the relevant dynamics including the nonlinear effects. Using the lowest number of frequencies possible (while still maintaining accuracy) is very important to minimize degrees of freedom and computation time for the optimization solver.\n",
+ "\n",
+ "The wave environment must be specified as a 2-dimensional `xarray.DataArray` containing the complex amplitude (m). \n",
+ "The two coordinates are the radial frequency ``omega`` (rad/s) and the direction ``wave_direction`` (rad). \n",
+ "The `wecopttool.waves` submodule contains functions for creating this `xarray.DataArray` for different types of wave environments. We will use the `wecopttool.waves.regular_wave` function to create the regular wave. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "amplitude = 0.0625 \n",
+ "wavefreq = 0.3\n",
+ "phase = 30\n",
+ "wavedir = 0\n",
+ "\n",
+ "f1 = wavefreq\n",
+ "nfreq = 10\n",
+ "\n",
+ "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Frequency and mesh check\n",
- "We will analyze 50 frequencies with a spacing of 0.05 Hz. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n",
+ "We already defined the fundamental frequency equal to the wave frequency and the number of frequencies to include all energetic harmonics. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n",
"\n",
"The `fb.minimal_computable_wavelength` parameter checks the mesh to determine the minimum wavelength that can be reliably computed using Capytaine. This warning is ignored here because the BEM results have been validated, but can be used as a guide for mesh refinement to ensure accurate BEM results."
]
@@ -131,8 +162,6 @@
"metadata": {},
"outputs": [],
"source": [
- "f1 = 0.05\n",
- "nfreq = 50\n",
"freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
"\n",
"min_computable_wavelength = fb.minimal_computable_wavelength\n",
@@ -252,32 +281,6 @@
"Note: We might receive a warning regarding negative linear damping values. Per default, WecOptTool ensures that the BEM data does not contain non-negative damping values. If you would like to correct the BEM solution manually to a minimum damping value you can specify `min_damping`. "
]
},
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Waves\n",
- "The wave environment must be specified as a 2-dimensional `xarray.DataArray` containing the complex amplitude (m). \n",
- "The two coordinates are the radial frequency ``omega`` (rad/s) and the direction ``wave_direction`` (rad). \n",
- "The `wecopttool.waves` submodule contains functions for creating this `xarray.DataArray` for different types of wave environments. \n",
- "\n",
- "In this case we will use a regular wave with a frequency of 0.3 Hz and an amplitude of 0.0625 m. \n",
- "We will use the `wecopttool.waves.regular_wave` function. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "amplitude = 0.0625 \n",
- "wavefreq = 0.3\n",
- "phase = 30\n",
- "wavedir = 0\n",
- "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
- ]
- },
{
"cell_type": "markdown",
"metadata": {},
@@ -820,8 +823,7 @@
"metadata": {},
"source": [
"### Results\n",
- "From a quick plot of the results, we see that the power absorption (where negative power is power absorbed by the device) generally improves for smaller values of `r2`.\n",
- "It is also clear that when the WEC is cylindrical (where `r2=0.88`), power absorption is reduced."
+ "From a quick plot of the results, we see that the power absorption (where negative power is power absorbed by the device) generally improves for larger values of `r2`."
]
},
{
diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
index 83df7e2d3..b97221d4d 100644
--- a/examples/tutorial_2_AquaHarmonics.ipynb
+++ b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -27,9 +27,6 @@
"import matplotlib.pyplot as plt\n",
"from matplotlib import cm\n",
"from scipy.optimize import brute\n",
- "import logging\n",
- "logging.basicConfig()\n",
- "logging.getLogger().setLevel(logging.DEBUG)\n",
"\n",
"import wecopttool as wot"
]
@@ -106,6 +103,33 @@
"fb.mass = np.atleast_2d(5e3) # [kg]"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Waves\n",
+ "A regular wave will allow us to get a good initial understanding of the optimal control trajectory.\n",
+ "Note that we need to choose an appropriate fundamental frequency, $f_1$, and number of frequencies, nfreq, to ensure that \n",
+ "the wave frequency and harmonic components are within the frequency array we use to calculate the hydrodynamic data."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "amplitude = 0.5\n",
+ "wavefreq = 0.24/2 \n",
+ "phase = 30\n",
+ "wavedir = 0\n",
+ "\n",
+ "f1 = wavefreq # [Hz]\n",
+ "nfreq = 10\n",
+ "\n",
+ "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -113,7 +137,7 @@
"#### Hydrodynamics\n",
"\n",
"Next we will run the boundary element model (BEM) [Capytaine](https://github.com/capytaine/capytaine) to obtain the hydrodynamic coefficients for the hull.\n",
- "Before running the BEM solver, we must specify a set of frequencies at which to perform the calculations.\n",
+ "Before running the BEM solver, we must create a set of frequencies at which to perform the calculations.\n",
"For `WecOptTool`, these must be a regularly spaced frequency array."
]
},
@@ -123,8 +147,6 @@
"metadata": {},
"outputs": [],
"source": [
- "f1 = 0.06 # [Hz]\n",
- "nfreq = 20\n",
"freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
"\n",
"bem_data = wot.run_bem(fb, freq, rho=rho, g=g)"
@@ -489,28 +511,6 @@
")"
]
},
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Waves\n",
- "A regular wave will allow us to get a good initial understanding of the optimal control trajectory.\n",
- "Note that we'll want to choose a wave frequency that is within the frequency array we used to calculate the hydrodynamic data."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "amplitude = 0.5\n",
- "wavefreq = 0.24/2 \n",
- "phase = 30\n",
- "wavedir = 0\n",
- "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
- ]
- },
{
"cell_type": "markdown",
"metadata": {},
diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
index 065ae977f..7bae07c5f 100644
--- a/examples/tutorial_3_LUPA.ipynb
+++ b/examples/tutorial_3_LUPA.ipynb
@@ -343,8 +343,18 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "### BEM\n",
- "With the LUPA geometry and physical properties fully defined, we can now run Capytaine to calculate the hydrodynamic and hydrostatic coefficients of the device, as done in previous tutorials. Capytaine can handle generalized modes and will calculate the coefficients for our 4 degrees of freedom. Because we are now using irregular waves, we need significantly more frequencies to capture the entire wave spectrum. The BEM coefficients have been pre-calculated and are saved in a file. To re-run the BEM, which takes about 1 hour, simply move or delete the existing `data/bem.nc` file."
+ "### Waves\n",
+ "Oregon State University has defined two sets of wave testing conditions for the LUPA: one corresponding to the PacWave South site and a scaling factor of 25, and one for the PacWave North site and a scaling factor of 20. \n",
+ "For each site/scale they provide four wave conditions to test at the Oregon State Large Wave Flume (LWF): the maximum 90th percentile, maximum percent annual energy, maximum occurrence, and minimum 10th percentile. \n",
+ "\n",
+ "In this tutorial we will use the PacWave South conditions scaled to the LWF and will design for the maximum occurrence wave. \n",
+ "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum. \n",
+ "\n",
+ "The irregular wave is created with multiple phase realizations. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes. For this tutorial, the number of realizations has been set to 2 to reduce the total run-time. \n",
+ "\n",
+ "$t_{total} = \\frac{nrealizations}{f1} $\n",
+ "\n",
+ "Because we are now using irregular waves, we need significantly more frequencies to capture the entire wave spectrum and WEC response. "
]
},
{
@@ -353,9 +363,57 @@
"metadata": {},
"outputs": [],
"source": [
- "# compute hydrodynamic coefficients\n",
+ "waves = {}\n",
+ "\n",
"f1 = 0.02\n",
"nfreq = 50\n",
+ "\n",
+ "# regular (for testing/setup)\n",
+ "amplitude = 0.1\n",
+ "wavefreq = 0.4\n",
+ "phase = 0\n",
+ "wavedir = 0\n",
+ "waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n",
+ "\n",
+ "nrealizations = 2\n",
+ "\n",
+ "# irregular wave cases from OSU\n",
+ "wave_cases = {\n",
+ " 'south_max_90': {'Hs': 0.21, 'Tp': 3.09}, \n",
+ " 'south_max_annual': {'Hs': 0.13, 'Tp': 2.35},\n",
+ " 'south_max_occurrence': {'Hs': 0.07, 'Tp': 1.90},\n",
+ " 'south_min_10': {'Hs': 0.04, 'Tp': 1.48}, \n",
+ " 'north_max_90': {'Hs': 0.25, 'Tp': 3.46}, \n",
+ " 'north_max_annual': {'Hs': 0.16, 'Tp': 2.63},\n",
+ " 'north_max_occurrence': {'Hs': 0.09, 'Tp': 2.13},\n",
+ " 'north_min_10': {'Hs': 0.05, 'Tp': 1.68}, \n",
+ "}\n",
+ "\n",
+ "def irregular_wave(hs, tp):\n",
+ " fp = 1/tp\n",
+ " spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
+ " efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
+ " return wot.waves.long_crested_wave(efth,nrealizations=nrealizations)\n",
+ "\n",
+ "for case, data in wave_cases.items():\n",
+ " waves[case] = irregular_wave(data['Hs'], data['Tp'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### BEM\n",
+ "With the LUPA geometry and physical properties fully defined, we can now run Capytaine to calculate the hydrodynamic and hydrostatic coefficients of the device, as done in previous tutorials. Capytaine can handle generalized modes and will calculate the coefficients for our 4 degrees of freedom. The BEM coefficients have been pre-calculated and are saved in a file. To re-run the BEM, which takes about 1 hour, simply move or delete the existing `data/bem.nc` file."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# compute hydrodynamic coefficients\n",
"freq = wot.frequency(f1, nfreq, False)\n",
"\n",
"# read BEM data file if it exists\n",
@@ -947,65 +1005,6 @@
")"
]
},
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Waves\n",
- "Oregon State University has defined two sets of wave testing conditions for the LUPA: one corresponding to the PacWave South site and a scaling factor of 25, and one for the PacWave North site and a scaling factor of 20. \n",
- "For each site/scale they provide four wave conditions to test at the Oregon State Large Wave Flume (LWF): the maximum 90th percentile, maximum percent annual energy, maximum occurrence, and minimum 10th percentile. \n",
- "\n",
- "In this tutorial we will use the PacWave South conditions scaled to the LWF and will design for the maximum occurrence wave. \n",
- "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum. \n",
- "\n",
- "The irregular wave is created with multiple phase realizations. For this tutorial, the number of realizations has been overwritten to 2 to reduce the total run-time. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes.\n",
- "\n",
- "$t_{total} = \\frac{nrealizations}{f1} $"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "waves = {}\n",
- "\n",
- "# regular (for testing/setup)\n",
- "amplitude = 0.1\n",
- "wavefreq = 0.4\n",
- "phase = 0\n",
- "wavedir = 0\n",
- "waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n",
- "\n",
- "# number of realizations to reach 20 minutes of total simulation time\n",
- "minutes_needed = 20\n",
- "nrealizations = minutes_needed*60*f1\n",
- "print(f'Number of realizations for a 20 minute total simulation time: {nrealizations}')\n",
- "nrealizations = 2 # overwrite nrealizations to reduce run-time\n",
- "\n",
- "# irregular wave cases from OSU\n",
- "wave_cases = {\n",
- " 'south_max_90': {'Hs': 0.21, 'Tp': 3.09}, \n",
- " 'south_max_annual': {'Hs': 0.13, 'Tp': 2.35},\n",
- " 'south_max_occurrence': {'Hs': 0.07, 'Tp': 1.90},\n",
- " 'south_min_10': {'Hs': 0.04, 'Tp': 1.48}, \n",
- " 'north_max_90': {'Hs': 0.25, 'Tp': 3.46}, \n",
- " 'north_max_annual': {'Hs': 0.16, 'Tp': 2.63},\n",
- " 'north_max_occurrence': {'Hs': 0.09, 'Tp': 2.13},\n",
- " 'north_min_10': {'Hs': 0.05, 'Tp': 1.68}, \n",
- "}\n",
- "\n",
- "def irregular_wave(hs, tp):\n",
- " fp = 1/tp\n",
- " spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
- " efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
- " return wot.waves.long_crested_wave(efth,nrealizations)\n",
- "\n",
- "for case, data in wave_cases.items():\n",
- " waves[case] = irregular_wave(data['Hs'], data['Tp'])"
- ]
- },
{
"cell_type": "markdown",
"metadata": {},
diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
index 4ab16ab2b..470e251ac 100644
--- a/examples/tutorial_4_Pioneer.ipynb
+++ b/examples/tutorial_4_Pioneer.ipynb
@@ -48,12 +48,12 @@
"## 1. Model setup\n",
"\n",
"### 1.1 Waves\n",
- "We start with the setting up the different waves we want to model, as this will inform what to select for our frequency range, which we need throughout the rest of the problem setup. \n",
+ "We start with setting up the different waves we want to model, as this will inform what to select for our frequency range, which we need throughout the rest of the problem setup. \n",
"We will consider two waves: a regular wave and an irregular wave, both with typical characteristics of the deployment site. The regular wave is roughly at 0.35 Hz, the known pitch resonance frequency of the buoy. The irregular wave has a peak period of 5 seconds, matching that of the deployment site. \n",
"\n",
- "The irregular wave is created with multiple phase realizations. For this tutorial, the number of realizations has been overwritten to 3 to reduce the total run-time. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations and average the resultant optimal power. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes.\n",
+ "Please refer to Tutorial 3 for information about selecting the number of realizations. For the purpose of the tutorial, the number of realizations has been set to 2 to reduce run-time.\n",
"\n",
- "$t_{total} = \\frac{nrealizations}{f1} $"
+ "The procedure for determining an appropriate frequency array for irregular wave conditions is detailed at the end of this tutorial."
]
},
{
@@ -62,8 +62,9 @@
"metadata": {},
"outputs": [],
"source": [
- "f1 = 0.025 # Hz\n",
- "nfreq = 75"
+ "fend = 1.875\n",
+ "nfreq = 150\n",
+ "f1 = fend/nfreq # Hz"
]
},
{
@@ -81,16 +82,12 @@
"Hs = 1.5\n",
"Tp = 5 \n",
"\n",
- "# number of realizations to reach 20 minutes of total simulation time\n",
- "minutes_needed = 20\n",
- "nrealizations = minutes_needed*60*f1\n",
- "print(f'Number of realizations for a 20 minute total simulation time: {nrealizations}')\n",
- "nrealizations = 3 # overwrite nrealizations to reduce run-time\n",
+ "nrealizations = 2\n",
"\n",
"fp = 1/Tp\n",
"spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, Hs)\n",
"efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
- "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations)"
+ "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations=nrealizations)"
]
},
{
@@ -854,7 +851,9 @@
"### 3.1 Solve\n",
"We will now run the same analysis for irregular waves. \n",
"\n",
- "An interesting result is that due to the narrow banded resonance of the flywheel, the controller attempts to make the excitation force monochromatic with the resonant frequency. To achieve this it uses significant reactive power (power by the PTO into the system). This is still worth it, resulting in a larger average electrical power output. "
+ "An interesting result is that due to the narrow banded resonance of the flywheel, the controller attempts to make the excitation force monochromatic with the resonant frequency. To achieve this it uses significant reactive power (power by the PTO into the system). This is still worth it, resulting in a larger average electrical power output. \n",
+ "\n",
+ "As noted previously, the optimization problem is solved for each wave realization. The optimal average power shown is the total average across the different realizations."
]
},
{
@@ -959,6 +958,119 @@
" axi.autoscale(axis='x', tight=True)"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 3.1 Notes on selection of frequency array and realizations\n",
+ "\n",
+ "In order to determine a suitable frequency range, it is best to complete an optimization with a large frequency range and examine the frequency domain results to see how large a range of frequencies is necessary to capture the relavent dynamics. Often, capturing nonlinear dynamics (such as due to constraints) requires a larger frequency range. In the case of the Pioneer WEC, it is clear that (likely due to the PTO torque constraint) the nonlinearities cause excitation at the odd harmonics of the excited frequency (which is slightly larger than the wave frequency itself for the Pioneer system). Specifically, the PTO force has significant peaks around the excited frequency and at the 3rd harmonic (with a small peak at the 5th harmonic as well) which leads to peaks in the frequency spectrum of electrical power at the 2nd, 4th, and 6th harmonics. The frequency values are normalized in the plots below to more clearly show the harmonics. In order to capture most of these nonlinearities, we chose a frequency range from 0 to 1.875 for this tutorial. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "filename = 'data/tutorial_4_freq_range.nc'\n",
+ "results = wot.read_netcdf(filename)\n",
+ "\n",
+ "excited_freq = abs(results['pos']).argmax()\n",
+ "normalized_freq = results['freq']/results['freq'][excited_freq]\n",
+ "\n",
+ "fig, ax = plt.subplots(5,1,figsize=(8, 12))\n",
+ "ax[0].stem(normalized_freq,abs(results['exc_force']))\n",
+ "ax[0].set_ylabel('Excitation Torque [N-m]')\n",
+ "ax[1].stem(normalized_freq,abs(results['pos']))\n",
+ "ax[1].set_ylabel('WEC Position [rad]')\n",
+ "ax[2].stem(normalized_freq,abs(results['vel']))\n",
+ "ax[2].set_ylabel('WEC Velocity [rad/s]')\n",
+ "ax[3].stem(normalized_freq,abs(results['pto_force']))\n",
+ "ax[3].set_ylabel('PTO Torque [N-m]')\n",
+ "ax[4].stem(normalized_freq,abs(results['power']))\n",
+ "ax[4].set_ylabel('Electrical Power [W]')\n",
+ "ax[4].set_xlabel('Frequency (normalized according to excited frequency)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "After determining a frequency range of interest, it is important to make sure enough frequencies are included in the array to generate accurate results. The Pioneer WEC has a very narrow banded response, which means it requires a large number of frequencies to model the response accurately. For any WecOptTool analysis in irregular wave conditions, it is best to complete a convergence study on the number of frequencies in the wave conditions of interest. The convergence study shown below varies the number of frequencies and includes enough realizations of each array to reach a 20 minute total simulation time (discussed further below). As shown below, an array of 150 frequencies is sufficient to model the system to within about 2% of the actual resultant mean power. Because the computation time increases with increasing number of frequencies, it is desirable to select the number of frequencies to minimize the computation time while meeting the intended accuracy. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "filename = 'data/tutorial_4_nfreqs.nc'\n",
+ "results = wot.read_netcdf(filename)\n",
+ "\n",
+ "power_percent_error = abs((results['power'] - results['power'][-1])/results['power'][-1])*100\n",
+ "time_per_realization = results['comp_time']/results['nrealizations']\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.plot(results['nfreqs'],power_percent_error)\n",
+ "plt.xlabel('Number of Frequencies in Array')\n",
+ "plt.ylabel('Absolute Mean Power Error [%]')\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.semilogy(results['nfreqs'],time_per_realization/60)\n",
+ "plt.xlabel('Number of Frequencies in Array')\n",
+ "plt.ylabel('Computation time per realization [min]')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, the number of phase realizations should be considered. Generally, it is recommended that enough realizations be used for the total simulation time to equal 20 minutes. \n",
+ "\n",
+ "$t_{total} = \\frac{nrealizations}{f1} $\n",
+ "\n",
+ "Although it usually leads to consistent results, this recommendation is somewhat arbitrary. To better understand the effect of the number of realizations on the overall result, it can be useful to complete a convergence study. The actual number of realizations needed depends on the problem itself (dynamics, constraints, etc.) and the desired precision of the result. As shown below, in the case of the Pioneer WEC, after about 15 realizations (20 minutes total simulation time) the result is within about 0.5% of the converged result which is deemed sufficient for this study."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "filename = 'data/tutorial_4_nrealizations.nc'\n",
+ "results = wot.read_netcdf(filename)\n",
+ "\n",
+ "rolling_mean = []\n",
+ "for ind in range(len(results['power'])):\n",
+ " rolling_mean.append(np.mean(results['power'][0:ind]))\n",
+ " \n",
+ "error_bar = rolling_mean[-1]*.005\n",
+ "\n",
+ "fig = plt.figure()\n",
+ "ax1 = fig.add_subplot(111)\n",
+ "ax2 = ax1.twiny()\n",
+ "ax1.plot(results['realization'],rolling_mean,label='Rolling Mean')\n",
+ "ax1.set_xlabel('Number of Realizations')\n",
+ "ax1.set_ylabel('Average Power [W]')\n",
+ "ax1.axhline(rolling_mean[-1]+error_bar,ls='--',label='0.5% error')\n",
+ "ax1.axhline(rolling_mean[-1]-error_bar,ls='--')\n",
+ "ax1.legend()\n",
+ "\n",
+ "tper_realization = (results['nfreq']/results['fend']).values # time per realization\n",
+ "\n",
+ "def tick_function(X):\n",
+ " V = tper_realization*X/60\n",
+ " return [\"%.1f\" % z for z in V]\n",
+ "\n",
+ "ax2.set_xticks(ax1.get_xticks())\n",
+ "ax2.set_xbound(ax1.get_xbound())\n",
+ "ax2.set_xticklabels(tick_function(ax1.get_xticks()))\n",
+ "ax2.set_xlabel('Total simulation time [min]')"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 18a74c250..6b2410451 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -326,10 +326,10 @@ def test_unstructured_controller_long_crested_wave(self,
res = wec.solve(waves=long_crested_wave,
obj_fun=pto.average_power,
nstate_opt=2*nfreq,
- x_wec_0=1e-1*np.ones(wec.nstate_wec),
- scale_x_wec=1e2,
- scale_x_opt=1e-2,
- scale_obj=1e-2,
+ x_wec_0=1e-3*np.ones(wec.nstate_wec),
+ scale_x_wec=1e1,
+ scale_x_opt=1e-3,
+ scale_obj=5e-2,
)
power_sol = -1*res[0]['fun']