diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index 45eed855b..60f9cc13d 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -5,14 +5,15 @@ Community contributions are welcomed! 🎊 * It is recommended that you create a *virtual environment*, e.g. using `conda`, `venv`, or similar. * If you want to build the documentation locally you will also need to install [pandoc](https://pandoc.org/installing.html) and [gifsicle](https://github.com/kohler/gifsicle). -* If you do not have Fortran compilers properly setup in your system, install `capytaine` and `wavesspectra` using `conda`. In this case you will need to have a `conda` virtual environment. This is recommended for *Windows* users since compiling `capytaine` on *Windows* requires [extra steps](https://github.com/capytaine/capytaine/issues/115). +* If you do not have Fortran compilers properly setup in your system, install `capytaine` and `wavespectra` using `conda`. In this case you will need to have a `conda` virtual environment. This is recommended for *Windows* users since compiling `capytaine` on *Windows* requires [extra steps](https://github.com/capytaine/capytaine/issues/115). +* Similarly, if you do not have the prerequisites listed [here](https://cyipopt.readthedocs.io/en/stable/install.html#from-source) to install `cyipopt` from source, install it using `conda` as well. * On a ZSH shell (*MacOS*) do `pip install -e .\[dev]` instead of `pip install -e .[dev]` in the instructions below (i.e., escape the opening square bracket). Using `conda` this looks like: ```bash conda create -n wecopttool python=3.11 conda activate wecopttool -conda install -c conda-forge capytaine wavespectra +conda install -c conda-forge capytaine wavespectra cyipopt git clone git@github.com:/WecOptTool.git cd WecOptTool pip install -e .[dev] @@ -42,7 +43,7 @@ Autograd does not support all NumPy and SciPy functionalities, see [autograd doc 1. Create a fork of WecOptTool 2. Create a branch for the specific issue 3. Add desired code modifications. For enhancements add to documentation. Add or modify a test. Make sure all tests pass and documentation builds. Follow style guide above. - 4. Do a pull request, and give admins edit access. [Link to any open issues](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue) and add relevant tags. Use a concise but descriptive PR title, as this will be part of the [release notes](https://github.com/sandialabs/WecOptTool/releases) for the next version. Start the PR title with an all caps label followed by a colon, e.g., "BUG FIX: ...", "NEW FEATURE: ...", "DOCUMENTATION: ...", etc. + 4. Do a pull request to the `dev` branch, and give admins edit access. [Link to any open issues](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue) and add relevant tags. Use a concise but descriptive PR title, as this will be part of the [release notes](https://github.com/sandialabs/WecOptTool/releases) for the next version. Start the PR title with an all caps label followed by a colon, e.g., "BUG FIX: ...", "NEW FEATURE: ...", "DOCUMENTATION: ...", etc. **Note: Pull requests should be made to the `dev` branch, not the `main` branch**. ## Tests There are a series of unit and integration tests defined in the `tests` directory. diff --git a/.github/RELEASING.md b/.github/RELEASING.md index 15d38664a..9148ec1ea 100644 --- a/.github/RELEASING.md +++ b/.github/RELEASING.md @@ -3,7 +3,8 @@ This section is for developers. Before a release make sure to: -* change [version number](https://semver.org/) in `pyproject.toml` +* change [version number](https://semver.org/) in `pyproject.toml` in the `dev` branch. +* Merge the `dev` branch into the `main` branch. **Note: the `dev` branch should only be merged into `main` when it is ready for a new release.** ## GitHub In the GitHub repository, click on *Releases*, click on *Draft new release*. diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 234d6f86b..a3fdb5d96 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -9,7 +9,7 @@ Concise description of the pull request with references to any [issues](https:// ## Checklist for PR - [ ] Authors read the [contribution guidelines](https://github.com/sandialabs/WecOptTool/blob/main/.github/CONTRIBUTING.md) -- [ ] The pull request is from an issue branch (not main) on *your* fork, to the [main branch in WecOptTool](https://github.com/sandialabs/WecOptTool). +- [ ] The pull request is from an issue branch (not main) on *your* fork, to the [dev branch in WecOptTool](https://github.com/sandialabs/WecOptTool/tree/dev). - [ ] The authors have given the admins edit access - [ ] All changes adhere to the style guide including PEP8, Docstrings, and Type Hints. - [ ] Modified the documentation if applicable diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index 28138c9da..f967e5dc4 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -4,9 +4,9 @@ name: "CodeQL" on: push: - branches: [ main ] + branches: [ main, dev ] pull_request: - branches: [ main ] + branches: [ main, dev ] schedule: - cron: '30 1 * * 0' diff --git a/.github/workflows/push.yml b/.github/workflows/push.yml index f8e3d0fee..e5e80be91 100644 --- a/.github/workflows/push.yml +++ b/.github/workflows/push.yml @@ -32,13 +32,15 @@ jobs: # - the content of pyproject.toml changes # - you manually change the value of the CACHE_NUMBER below # Else the existing cache is used. - - name: Setup Mambaforge - uses: conda-incubator/setup-miniconda@v2 + - name: Setup Miniforge + uses: conda-incubator/setup-miniconda@v3 with: - miniforge-variant: Mambaforge - miniforge-version: latest - activate-environment: test-env - use-mamba: true + miniforge-version: latest + mamba-version: "*" + python-version: '3.11' + channels: conda-forge + channel-priority: true + activate-environment: test-env # save the date to include in the cache key - name: Get Date @@ -57,6 +59,7 @@ jobs: echo " - python=${{ matrix.python-version }}" >> environment.yml echo " - capytaine=1.5" >> environment.yml echo " - wavespectra" >> environment.yml + echo " - cyipopt" >> environment.yml cat environment.yml # use the cache if it exists @@ -127,6 +130,16 @@ jobs: - uses: actions/setup-python@v4 with: python-version: '3.11' # CHANGE: Python version + + - name: Setup Miniforge + uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + mamba-version: "*" + python-version: '3.11' + channels: conda-forge + channel-priority: true + activate-environment: test-env - name: Install dependencies run: sudo apt-get install libglu1 pandoc gifsicle @@ -134,6 +147,7 @@ jobs: - name: Install WecOptTool for documentation shell: bash -l {0} run: | + mamba install cyipopt python3 -m pip install --upgrade pip pip3 install wheel coveralls pip3 install .[dev,geometry] diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 72a979038..8a7bd0f50 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -49,12 +49,23 @@ jobs: with: python-version: '3.11' + - name: Setup Miniforge + uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + mamba-version: "*" + python-version: '3.11' + channels: conda-forge + channel-priority: true + activate-environment: test-env + - name: Install dependencies run: sudo apt-get install libglu1 pandoc gifsicle - name: Install WecOptTool for documentation shell: bash -l {0} run: | + mamba install cyipopt python3 -m pip install --upgrade pip pip3 install wheel coveralls pip3 install .[dev,geometry] 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/docs/source/wecopttool_refs.bib b/docs/source/wecopttool_refs.bib index fe4cf2645..21e553e2d 100644 --- a/docs/source/wecopttool_refs.bib +++ b/docs/source/wecopttool_refs.bib @@ -7,6 +7,7 @@ @inproceedings{Gaebele:2023wf date-added = {2023-01-12 14:02:25 -0700}, date-modified = {2023-10-03 08:41:16 -0600}, month = {June}, + doi = {10.1115/OMAE2023-103899}, series = {International Conference on Ocean, Offshore and Arctic Engineering}, title = {INCORPORATING EMPIRICAL NONLINEAR EFFICIENCY INTO CONTROL CO-OPTIMIZATION OF A REAL WORLD HEAVING POINT ABSORBER USING {WECOPTTOOL}}, year = {2023}} 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..b69497292 100644 --- a/examples/tutorial_1_WaveBot.ipynb +++ b/examples/tutorial_1_WaveBot.ipynb @@ -6,11 +6,11 @@ "source": [ "# Tutorial 1 - WaveBot\n", "The goal of this tutorial is to familiarize new users with how to set up and run optimization problems using WecOptTool. \n", - "It uses a one-body WEC, the WaveBot, in one degree of freedom in regular waves. \n", + "This tutorial uses a one-body WEC, the WaveBot, in one degree of freedom in regular waves. \n", "\n", "![WaveBot Photo](https://live.staticflickr.com/65535/51855905347_de87ccaaba_z.jpg)\n", "\n", - "At the end of this tutorial the user will perform control co-design of the WEC's geometry and a corresponding optimal controller to maximize electrical power. \n", + "At the end of this tutorial the user will perform control co-design of the WEC geometry and a corresponding optimal controller to maximize electrical power. \n", "We build up to this problem in three parts of successive complexity:\n", "\n", "1. [Optimal control for maximum mechanical power](#1.-Optimal-control-for-maximum-mechanical-power)\n", @@ -19,8 +19,8 @@ "\n", "We will start by loading the necessary modules: \n", "\n", - "* Import Autograd (wrapper on NumPy, required) for automatic differentiation\n", - "* Import other packages we will use in this tutorial \n", + "* Import Autograd (wrapper on NumPy, required) for [automatic differentiation](https://sandialabs.github.io/WecOptTool/theory.html#automatic-differentiation)\n", + "* Import other packages we will use in this tutorial\n", "* Import WecOptTool " ] }, @@ -35,7 +35,10 @@ "import matplotlib.pyplot as plt\n", "from scipy.optimize import brute\n", "\n", - "import wecopttool as wot" + "import wecopttool as wot\n", + "\n", + "## set colorblind-friendly colormap for plots\n", + "plt.style.use('tableau-colorblind10')" ] }, { @@ -43,42 +46,96 @@ "metadata": {}, "source": [ "## 1. Optimal control for maximum mechanical power\n", - "This example illustrates how to set up, run, and analyze a basic optimization problem within WecOptTool.\n", + "This example illustrates how to set up, run, and analyze a basic optimization problem using WecOptTool.\n", "\n", - "The objective of this example is to **find the optimal PTO force time-series** that produces the most mechanical power subject to the WEC dynamics and a maximum force the PTO can exert.\n", + "The objective of this example is to **find the optimal PTO force time series** that produces the most mechanical power subject to the WEC dynamics and a maximum force the PTO can exert.\n", "\n", "WecOptTool requires the following to be defined to successfully run its optimization routines:\n", - "- The WEC object, including all of its properties and constraints\n", - "- The wave condition\n", - "- The objective function\n", + "\n", + "* The wave environment\n", + "* The WEC object, including all of its properties and constraints\n", + "* The objective function\n", + "\n", + "The graphic below shows all the requirements for this first part of the tutorial: the wave environment on the left, the WEC object(with its relevant subcomponents) in the blue box in the middle, and the objective function (mechanical power) on the right.\n", "\n", "
\n", "\n", - "
\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Waves and WEC geometry\n", + "The pseudo-spectral solution used in WecOptTool requires a specified set of frequencies to analyze.\n", + "To model the WEC accurately, we need to ensure the set of selected frequencies captures the full hydrodynamic response range of the WEC.\n", + "\n", + "Therefore, we start our WecOptTool model by defining this set of frequencies, the wave environment, and the mesh of the WEC geometry.\n", + "This ensures the minimum wavelength in the selected frequencies is larger than the minimum wavelength that can be analyzed with the mesh using the Boundary Element Method (BEM), which we will do [later](#bem-solution).\n", + "This is a good initial check to make sure the hydrodynamic properties results calculated with the BEM will be accurate.\n", + "\n", + "#### Frequency selection\n", + "To start with a simple case in this tutorial, we are interested in modeling a regular wave at 0.3 Hz.\n", + "In regular waves, the linear WEC system response will always occur at the wave frequency, with the nonlinear system response also occurring at the corresponding odd harmonics (3rd, 5th, etc.) of the wave frequency.\n", + "Thus, we can set the fundamental frequency in our set, $f_1$, equal to the wave frequency.\n", + "Prior WaveBot analysis has found that the nonlinear dynamics are sufficiently captured with the 3rd, 5th, 7th, and 9th harmonics.\n", + "WecOptTool assumes the spacing between frequencies equals the magnitude of $f_1$. Therefore, we will model 10 frequencies to encompass these harmonics.\n", "\n", - "The graphic shows all the requirements for this first part of the tutorial: from the wave on the left, to the objective (mechanical power) on the right.\n", - "The WEC object, with all it's components, is illustrated in the middle. The components inside the blue box are the WEC properties that are actually passed on to the optimizer.\n", - "In short, the WEC's hydrodynamic properties are modelled by\n", - "1. Defining the WEC's geometry\n", - "2. Meshing the geometry\n", - "3. Obtaining the WEC's BEM cofficients based on the mesh\n", - "4. Determening the WEC's intrinsic impedance model based on the BEM coefficients\n", + "These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem.\n", + "See the [Theory](https://sandialabs.github.io/WecOptTool/theory.html) section of the WecOptTool documentation for more details on the pseudo-spectral problem formulation.\n", "\n", - "For this first part of the tutorial, the heave-only, WEC-PTO kinematics are trivial (Unity) and the PTO is assumed to be lossless." + "It is important to use the lowest number of frequencies possible while still maintaining accuracy in order to minimize degrees of freedom and computation time for the optimization solver." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "### WEC object\n", - "In this section we will create the `WEC` object, which contains all the information about the WEC and its dynamics. This constitutes the vast majority of the setup required to run WecOptTool.\n", + "wavefreq = 0.3 # Hz\n", + "f1 = wavefreq\n", + "nfreq = 10\n", "\n", - "Our `WEC` object requires information about the mesh, degrees of freedom, mass and hydrostatic properties, linear hydrodynamic coefficients (from a BEM solution), any additional dynamic forces (e.g. PTO force, mooring, non-linear hydrodynamics), and constraints (e.g. maximum PTO extension). \n", - "In this case, the only additional force will be the PTO force and the only constraint will be a maximum PTO force of $2,000 N$.\n", + "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Waves\n", + "WecOptTool is configured to have the wave environment specified as a 2-dimensional `xarray.DataArray` containing the complex wave amplitudes (in meters), wave phases (in degrees), and directions (in degrees).\n", + "We will use an amplitude of 0.0625 meters, a phase of 30 degrees, and direction 0 for this tutorial.\n", + "The two coordinates are the radial frequency ``omega`` (in rad/s) and the direction ``wave_direction`` (in radians).\n", + "\n", + "The `xarray` package has a pretty steep learning curve, so WecOptTool includes the `waves` module containing more intuitive functions that create `xarray.DataArray` objects for different types of wave environments.\n", + "See the `xarray` documentation [here](https://docs.xarray.dev/en/latest/index.html) if you are interested in learning more about the package.\n", + "In this case, we will use the `wecopttool.waves.regular_wave` function with our wave frequency, amplitude, phase, and direction:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "amplitude = 0.0625 # m\n", + "phase = 30 # degrees\n", + "wavedir = 0 # degrees\n", "\n", - "#### Mesh\n", - "First, we will create a surface mesh for the hull and store it using the `FloatingBody` object from Capytaine. The WaveBot mesh is pre-defined in the `wecopttool.geom` module, so we will call it directly from there. We will only model the heave degree of freedom in this case. Note that the Capytaine `from_meshio` method can also import from other file types (STL, VTK, MSH, etc.)" + "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### WEC geometry mesh\n", + "Now we will create a surface mesh for the WEC hull and store it using the `FloatingBody` object from Capytaine.\n", + "The WaveBot mesh is pre-defined in the `wecopttool.geom` module, so we will call it directly from there.\n", + "Note that the Capytaine `from_meshio` method can also import from other file types (STL, VTK, MSH, etc.), click [here](https://pypi.org/project/meshio/) for the full list of compatible mesh file types." ] }, { @@ -90,19 +147,17 @@ "wb = wot.geom.WaveBot() # use standard dimensions\n", "mesh_size_factor = 0.5 # 1.0 for default, smaller to refine mesh\n", "mesh = wb.mesh(mesh_size_factor)\n", - "fb = cpy.FloatingBody.from_meshio(mesh, name=\"WaveBot\")\n", - "fb.add_translation_dof(name=\"Heave\")\n", - "ndof = fb.nb_dofs" + "fb = cpy.FloatingBody.from_meshio(mesh, name=\"WaveBot\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "At this point we can visualize the mesh for inspection.\n", + "We can also visualize the mesh for inspection.\n", "Capytaine has built-in methods for visualizing meshes (`fb.show`, and `fb.show_matplotlib`). \n", - "When running outside a Notebook, these are interactive. \n", - "The included WaveBot example also has a method for plotting the cross-section of the device. " + "These methods are interactive when used outside a Jupyter notebook. \n", + "The included WaveBot example also has a method for plotting the cross-section of the device, `plot_cross_section`. " ] }, { @@ -112,17 +167,19 @@ "outputs": [], "source": [ "fb.show_matplotlib()\n", - "_ = wb.plot_cross_section(show=True) # specific to WaveBot" + "_ = wb.plot_cross_section(show=True) # specific to WaveBot" ] }, { "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", - "\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." + "#### Minimum wavelength check\n", + "With the frequency vector, wave environment, and geometry mesh all defined, we can now check to make sure they are all suitable to accurately simulate the WEC dynamics.\n", + "The `fb.minimal_computable_wavelength` method checks the mesh to determine the minimum wavelength that can be reliably computed using Capytaine.\n", + "We compare this value to the minimum wavelength our frequency vector will compute; we want this number to be larger than Capytaine's minimum wavelength.\n", + "A warning is printed if this is not the case.\n", + "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,10 +188,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", "g = 9.81\n", "min_period = 1/(f1*nfreq)\n", @@ -149,10 +202,83 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### BEM\n", - "With our Capytaine floating body created, we can now run the Boundary Element Method solver in Capytaine to get the hydrostatic and hydrodynamic coefficients of our WEC object. This is wrapped into the `wecopttool.run_bem` function.\n", + "### WEC object\n", + "The `WEC` object in WecOptTool contains all the information about the WEC and its properties and dynamics.\n", + "This constitutes the vast majority of the setup required to run a WecOptTool optimization.\n", + "The WEC object handles three categoreis of data that are passed on to the optimizer:\n", + "\n", + "1. The intrinsic impedance of the WEC\n", + "2. The WEC and PTO kinematic\n", + "3. The PTO dynamics\n", + "\n", + "In order for the `WEC` object to be able to compute these data, we must provide information when we declare the object in the code.\n", + "The required information includes:\n", + "\n", + "* [Degrees of freedom to consider](#degrees-of-freedom)\n", + "* [Linear hydrodynamic coefficients](#bem-solution)\n", + "* Any additional loads (e.g. PTO force, mooring, non-linear hydrodynamics)\n", + "* Device constraints\n", + "\n", + "Again, we will keep things simple to start.\n", + "We will only consider heave motion and assume a lossless PTO, and the WaveBot has trivial WEC-PTO kinematics (unity).\n", + "We will apply one additional force (the PTO force on the WEC) and one constraint (the maximum PTO force), which [we define below using the PTO class](#pto)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Degrees of freedom\n", + "The degrees of freedom can be added directly to the `FloatingBody` using the `add_translation_dof` and `add_rotation_dof` methods. The axis of rotation is automatically defined if the method is supplied a name of `\"Surge\"`, `\"Sway\"`, `\"Heave\"`, `\"Roll\"`, `\"Pitch\"`, or `\"Yaw\"`. We will only model the heave degree of freedom in this case." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fb.add_translation_dof(name=\"Heave\")\n", + "ndof = fb.nb_dofs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### BEM solution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can use the `FloatingBody` and frequency vector created earlier to run the Boundary Element Method solver in Capytaine to get the hydrostatic and hydrodynamic coefficients of our WEC object.\n", + "This is wrapped into the `wecopttool.run_bem` function.\n", + "\n", + "If you would like to save the BEM data to a NetCDF file for future use, uncomment the second line of the cell below to use the `wot.write_netcdf` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bem_data = wot.run_bem(fb, freq)\n", + "# wot.write_netcdf('bem_data.nc', bem_data) # saves BEM data to file" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we use the utilities class `wecopttool.utilities`, which has functions that help you analyze and design WECs, but are not needed for the core function of WecOptTool.\n", + "For example, we calculate the WEC's intrinsic impedance with its hydrodynamic coefficients and inherent inertial properties. We make use of `wecopttool.add_linear_friction()` to convert the `bem_data` into a dataset which contains all zero friction data, because `wecopttool.check_radiation_damping()` and `wecopttool.hydrodynamic_impedance()` excpect a data variable called 'friction'.\n", "\n", - "If you would like to save our BEM data to a NetCDF file for future use, see the `wecopttool.write_netcdf` function." + "The intrinsic impedance tells us how a hydrodynamic body responds to different excitation frequencies. For oscillating systems we are intersted in both, the amplitude of the resulting velocity as well as the phase between velocity and excitation force. Bode plots are useful tool to visualize the frequency response function.\n", + "\n", + "The natural frequency reveals itself as a trough in the intrinsic impedance for restoring degrees of freedom (heave, roll, pitch)." ] }, { @@ -161,7 +287,13 @@ "metadata": {}, "outputs": [], "source": [ - "bem_data = wot.run_bem(fb, freq)" + "hd = wot.add_linear_friction(bem_data, friction = None) \n", + "#we're not actually adding friction, but need the datavariables in hd \n", + "hd = wot.check_radiation_damping(hd)\n", + "\n", + "intrinsic_impedance = wot.hydrodynamic_impedance(hd)\n", + "fig, axes = wot.utilities.plot_bode_impedance(intrinsic_impedance,\n", + " 'WaveBot Intrinsic Impedance')" ] }, { @@ -169,14 +301,20 @@ "metadata": {}, "source": [ "#### PTO\n", - "WecOptTool includes the `PTO` class to encompass all properties of the power take-off system of the WEC. Data wrapped into our `PTO` class will be used to help define our `WEC` object and optimization problem later.\n", + "WecOptTool includes the `PTO` class to encompass all properties of the power take-off system of the WEC.\n", + "Data wrapped into our `PTO` class will be used to provide the additional force and constraint for our `WEC` object definition, as well as our optimization problem later.\n", "\n", "To create an instance of the `PTO` class, we need:\n", - "- The kinematics matrix, which converts from the WEC degrees of freedom to the PTO degrees of freedom. The PTO extracts power directly from the WEC's heave in this case, so the kinematics matrix is simply the $1 \\times 1$ identity matrix.\n", - "- The definition of the PTO controller. The `wecopttool.pto` submodule includes P, PI, and PID controller functions that can be provided to the `PTO` class and return the PTO force. However, we will be using an unstructured controller in this case, so we will set `None` for the controller.\n", - "- Any PTO impedance. We're only interested in mechanical power for this first problem, so we will leave this empty for now\n", - "- The non-linear power conversion loss (assumed 0% if `None`)\n", - "- The PTO system name, if desired" + "\n", + "* The kinematics matrix, which converts from the WEC degrees of freedom to the PTO degrees of freedom.\n", + "In this case, the PTO extracts power directly from the heaving motion of the WEC, so the kinematics matrix is simply a $1 \\times 1$ identity matrix.\n", + "* The definition of the PTO controller.\n", + "The `wecopttool.pto` submodule includes P, PI, and PID controller functions that can be provided to the `PTO` class and return the PTO force.\n", + "However, we will be using an unstructured controller in this case, which is the default when `controller=None`.\n", + "* Any PTO impedance.\n", + "We are only interested in mechanical power for this first problem, so we will leave this empty for now.\n", + "* The non-linear power conversion loss (assumed 0% if `None`)\n", + "* The PTO system name, if desired" ] }, { @@ -197,7 +335,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's define the PTO forcing on the WEC and the PTO constraints. For our optimization problem, the constraints must be in the correct format for `scipy.optimize.minimize()`. We will enforce the constraint at 4 times more points than the dynamics (see Theory for why this is helpful for the pseudo-spectral problem)." + "Now we define the PTO forcing on the WEC and the PTO constraints.\n", + "We set the maximum PTO force as 750 Newtons.\n", + "This value is chosen somewhat arbitrary to highlight different effects throughout this tutorial.\n", + "In practice, the user would need to identify their limiting component in the PTO and then compute a suitable value.\n", + "For example, if the generator drive has a maximal current $I_{max}$ of 10 A, with a generator torque constant $K_t$ of 8 Nm/A and a gear ratio $N$ of 9 rad/m, this results in a max PTO linear force of $F = N K_{t} I_{max} = 9$ rad/m $\\times 8$ Nm $\\times 10$ A $= 720$ N.\n", + "\n", + "We will enforce the constraint at 4 times more points than the dynamics; see [Theory]((https://sandialabs.github.io/WecOptTool/theory.html)) for why this is helpful for the pseudo-spectral problem.\n", + "The constraints must be in the correct format for `cyipopt.minimize_ipopt`, which is the solver WecOptTool uses. See the documentation for the solver [here](https://cyipopt.readthedocs.io/en/stable/reference.html#cyipopt.minimize_ipopt) and note the matching syntax of `ineq_cons` below. Note that this is the same syntax as used in old versions of WecOptTool for `scipy.optimize.minimize`." ] }, { @@ -210,10 +355,10 @@ "f_add = {'PTO': pto.force_on_wec}\n", "\n", "# Constraint\n", - "f_max = 2000.0\n", + "f_max = 750.0\n", "nsubsteps = 4\n", "\n", - "def const_f_pto(wec, x_wec, x_opt, waves): # Format for scipy.optimize.minimize\n", + "def const_f_pto(wec, x_wec, x_opt, waves): # Format for cyipopt.minimize_ipopt\n", " f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n", " return f_max - np.abs(f.flatten())\n", "\n", @@ -228,7 +373,9 @@ "metadata": {}, "source": [ "#### `WEC` creation\n", - "We are now ready to create the `WEC` object itself! Since we ran our BEM already, we can define the object using the `wecopttool.WEC.from_bem` function. If we saved our BEM data to a NetCDF file, we can also provide the path to that file instead of specifying the BEM `Dataset` directly." + "We are now ready to create the `WEC` object itself!\n", + "Since we ran the BEM already, we can define the object using the `wecopttool.WEC.from_bem` function.\n", + "If we saved BEM data to a NetCDF file, we can also provide the path to that file instead of specifying the BEM `Dataset` directly." ] }, { @@ -249,20 +396,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "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`. " + "_Note: You may receive a warning regarding negative linear damping values._\n", + "_By default, WecOptTool ensures that the BEM data does not contain non-negative damping values._\n", + "_If you would like to correct the BEM solution to a different minimum damping value besides zero, 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. " + "### Objective function\n", + "The objective function is the quantity we want to optimize—in this case, the average mechanical power.\n", + "The function to compute this can be taken directly from the `PTO` object we created:" ] }, { @@ -271,21 +416,17 @@ "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)" + "obj_fun = pto.mechanical_average_power" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Objective function\n", - "The objective function is the quantity (scalar) we want to optimize—in this case, the average mechanical power. The objective function is itself a function of the optimization state, the size of which we need to properly define our call to `scipy.optimize.minimize()`. The average mechanical power can be taken directly from the `PTO` object we created.\n", - "\n", - "One technical quirk here: `nstate_opt` is one smaller than would be expected for a state space representing the mean (DC) component and the real and imaginary Fourier coefficients. This is because WecOptTool excludes the imaginary Fourier component of the highest frequency (the 2-point wave). Since the 2-point wave is sampled at multiples of $\\pi$, the imaginary component is evaluated as $sin(n\\pi); n = 0, 1, 2, ..., n_{freq}$, which is always zero. Excluding this component speeds up the optimization as the state space is reduced by one." + "The objective function is itself a function of the optimization state `x_opt`.\n", + "The length of `x_opt`, `nstate_opt`, needs to be properly defined to successfully call `cyipopt.minimize_ipopt`.\n", + "In other words, it should match the quantities we are interested in optimizing.\n", + "In this case, we are optimizing for the control trajectories of an unstructured PTO, which can be represented in the Fourier domain by the DC (zero frequency) component, then the real and imaginary components for each frequency." ] }, { @@ -294,23 +435,43 @@ "metadata": {}, "outputs": [], "source": [ - "obj_fun = pto.mechanical_average_power\n", "nstate_opt = 2*nfreq" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One technical quirk here: `nstate_opt` is one smaller than would be expected for a state space representing the mean (DC) component and the real and imaginary Fourier coefficients.\n", + "This is because WecOptTool excludes the imaginary Fourier component of the highest frequency (the 2-point wave).\n", + "Since the 2-point wave is sampled at multiples of $\\pi$, the imaginary component is evaluated as $sin(n\\pi); n = 0, 1, 2, ..., n_{freq}$, which is always zero." + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve\n", - "We are now ready to solve the problem. WecOptTool uses `scipy.optimize.minimize` as its optimization driver, which is wrapped into `wecopttool.WEC.solve` for ease of use.\n", + "We are now ready to solve the problem.\n", + "WecOptTool uses `cyipopt.minimize_ipopt` as its optimization driver, which is wrapped into `wecopttool.WEC.solve` for ease of use.\n", "\n", - "Note that the only required inputs for defining and solving the problem are: (1) the waves, (2) the objective function, and (3) the size of the optimization state. Optional inputs can be provided to control the optimization execution if desired, which we do here to change the default iteration maximum and tolerance. See `scipy.optimize.minimize` docs for more details.\n", + "The only required inputs for defining and solving the problem are:\n", "\n", - "To help the optimization we will scale the problem before solving it (see Documentation). WecOptTool allows you to scale the WEC dynamics state, your optimization state (in this case the Fourier coefficients for the PTO force), and the objective function separately. See the `wecopttool.WEC.solve()` function for more information.\n", + "1. The wave environment\n", + "2. The objective function\n", + "3. The size of the optimization state (`nstate_opt`)\n", "\n", + "Optional inputs can be provided to control the optimization execution if desired, which we do here to change the default iteration maximum and tolerance.\n", + "See the `cyipopt.minimize_ipopt` documentation [here](https://cyipopt.readthedocs.io/en/stable/reference.html#cyipopt.minimize_ipopt) for more details.\n", "\n", - "Pay attention to the `Exit mode`: an exit mode of $0$ indicates a successful solution. For an easy problem (linear, single Dof, unconstrained, etc.) your iterations shouldn't need to exceed 100. If they do, try adjusting the scales by orders of magnitude, one at a time." + "To help the problem converge faster, we will scale the problem before solving it (see the [Scaling](https://sandialabs.github.io/WecOptTool/theory.html#scaling) section of the Theory documentation).\n", + "WecOptTool allows you to scale the WEC dynamics state, the optimization state, and the objective function separately.\n", + "See the `wecopttool.WEC.solve` documentation [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.WEC.solve.html#wecopttool-wec-solve).\n", + "\n", + "\n", + "Pay attention to the `Exit mode`: an exit mode of `0` indicates a successful solution.\n", + "A simple problem (linear, single degree of freedom, unconstrained, etc.) should converge in well under 100 iterations.\n", + "If you exceed this, try adjusting the scales by orders of magnitude, one at a time." ] }, { @@ -343,7 +504,10 @@ "metadata": {}, "source": [ "### Analyzing results\n", - "We will use two post-processing functions to obtain frequency- and time-domain results for the WEC and PTO responses. The pseudospectral method gives continuous in time results. To get smoother looking plots, we specify the number of subpoints betweeen co-location points. In this case we will use 5. " + "We will use two post-processing functions to obtain frequency and time domain results for the WEC and PTO responses.\n", + "The `wec.post_process` and `pto.post_process` functions both post-process the results for each wave phase realization.\n", + "The pseudo-spectral method gives continuous in time results.\n", + "To get smoother looking plots, we specify five subpoints betweeen co-location points." ] }, { @@ -353,15 +517,16 @@ "outputs": [], "source": [ "nsubsteps = 5\n", - "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves.sel(realization=0), nsubsteps=nsubsteps)\n", - "wec_fdom, wec_tdom = wec.post_process(results[0], waves.sel(realization=0), nsubsteps=nsubsteps)" + "pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps)\n", + "wec_fdom, wec_tdom = wec.post_process(wec, results, waves, nsubsteps=nsubsteps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The `pto.post_process` function returns `xarray.Dataset`s, which have built-in integration with PyPlot for smart plotting that automagically sets titles and formatting. We will plot the mechanical power (`mech_power`), position (`pos`), and the PTO force (`force`)." + "The `pto.post_process` function returns a list of `xarray.Dataset` objects, each element of which has built-in integration with PyPlot for smart plotting with automatic titles and formatting.\n", + "We will plot the mechanical power (`mech_power`), position (`pos`), and the PTO force (`force`)." ] }, { @@ -371,14 +536,15 @@ "outputs": [], "source": [ "plt.figure()\n", - "pto_tdom['power'].loc['mech',:,:].plot()" + "pto_tdom[0]['power'].loc['mech',:,:].plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We could similarly plot any time- or frequency-domain repsonse of the WEC or PTO. For instance, here is the PTO heave motion." + "We could similarly plot any time or frequency domain response of the WEC or PTO by calling the specific type of response (position, velocity, force, etc.).\n", + "For example, to plot the WEC heave position and PTO force:" ] }, { @@ -388,7 +554,18 @@ "outputs": [], "source": [ "plt.figure()\n", - "wec_tdom['pos'].plot()" + "wec_tdom[0]['pos'].plot()\n", + "\n", + "plt.figure()\n", + "pto_tdom[0]['force'].plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that there are other dynamic responses available in the post-processed WEC and PTO variables (`wec_tdom`, `pto_tdom`, `wec_fdom`, `pto_fdom`).\n", + "For example, the time domain PTO variable contains the following response:" ] }, { @@ -397,15 +574,14 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure()\n", - "pto_tdom['force'].plot()" + "pto_tdom" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Note that there are other dynamic responses available in the post-processed WEC and PTO variables (`wec_tdom`, `pto_tdom`, `wec_fdom`, `pto_fdom`). For example, the time domain PTO variable contains the following response:" + "Lastly, we will visualize the average power at different stages and how the power flows through the system.\n" ] }, { @@ -414,7 +590,31 @@ "metadata": {}, "outputs": [], "source": [ - "pto_tdom" + "p_flows = wot.utilities.calculate_power_flows(wec,\n", + " pto, \n", + " results, \n", + " waves, \n", + " intrinsic_impedance)\n", + "wot.utilities.plot_power_flow(p_flows)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "On the very left of the power flow Sankey diagram we start with the *Optimal Excitation*, which is only determined by the incident wave and the hydrodynamic damping and friction. \n", + "In order for the actual *Excited* power to equal the *Optimal Excitation* the absorbed *Mechanical* power would need to equal the *Radiated* power (power that is radiated away from the WEC). In this case the *Unused Potential* would be zero.\n", + "In other words, you can never convert more than 50% of the incident wave energy into kinetic energy. For more information on this we refer to Johannes Falnes Book - Ocean Waves and Oscillating System, specifically Chapter 6.\n", + "\n", + "However, the optimal 50% absorption is an overly optimistic goal considering that real world systems are likely constrained in their oscillation velocity amplitude and phase, due to limitations in stroke and applicable force.\n", + "\n", + "\n", + "The PTO force constraint used in this optimization also stopped us from absorbing the maximal potential wave energy. It is notable that we absorbed approximately 3/4 of the maximal absorbable power (*Mechanical* / (1/2 *Optimal Excitation*)), with relatively little *Radiated* power (about 1/2 of the absorbed power). To also absorb the last 1/4 of the potential wave power, we would need to increase the *Radiated* power three times!!\n", + "\n", + "\n", + "A more important question than \"How to achieve optimal absorption?\" is \"How do we optimize usable output power?\", i.e. *Electrical* power. For this optimization case we used a lossless PTO that has no impedance in itself. Therefore, the *Electrical* power equals the absorbed *Mechanical* power.\n", + "We'll show in the following sections that this a poor assumption and that the power flow looks fundamentally different when taking the PTO dynamics into account!!\n" ] }, { @@ -423,19 +623,21 @@ "source": [ "## 2. Optimal control for maximum electrical power\n", "\n", - "The rest of this tutorial will focus on optimizing for electrical power (new objective function) rather than mechanical, as this is a form of power that is usable and transportable.\n", + "The rest of this tutorial will focus on a different objective function: we will now optimize for electrical power rather than mechanical, as this is a form of power that is usable and transportable.\n", "\n", - "Since we're still dealing with the same WaveBot as in part 1, we can reuse the BEM and wave data from before. Look back at part 1 if you need a refresher on how to create these data.\n", + "Since we are still dealing with the same WaveBot as in Part 1, we can reuse the BEM, WEC-PTO kinematics, and wave data from before.\n", "\n", "
\n", "\n", "
\n", - "The WEC-PTO kinematics remain the same as well (unity). The major difference now is that we consider the dynamics of PTO, since they impact the electrical power and we shall not assume a lossless PTO.\n", "\n", - "We will express the PTO's dynamics in form of a 2-port impedance model, to incoporate the dynamics of the drive-train and the dynamics of the generator.\n", - "The additional mechanical energy storage through the drive-train is modelled using Newton's second law and we assume a linear generator using a power-invariant park transform.\n", + "The biggest difference when considering electrical power is the addition of the PTO dynamics.\n", + "In other words, we no longer assume a lossless PTO.\n", + "We will express the PTO dynamics in the form of a [2-port impedance model](https://en.wikipedia.org/wiki/Impedance_parameters) to incorporate the dynamics of the drivetrain and generator.\n", + "The additional mechanical energy storage through the drivetrain is modeled using Newton's second law, and we assume a linear generator using a power-invariant park transform.\n", "\n", - "The PTO impedance matrix components are then obtained under open-circuit conditions, i.e. no load current or no WEC velocity, respectively." + "The PTO impedance matrix components are then obtained under open-circuit conditions: no load current ($i_{pto}$ in figure above) or no WEC velocity ($vel_{pto}$ above), respectively.\n", + "See Section II of Michelén et al. [here](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=10114969) for more details about how the code below is derived." ] }, { @@ -474,7 +676,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, we will create a new `PTO` object with this impedance matrix. We will also update the definitions of our PTO constraint and additional dynamic forcing function to use the new object. We will set our PTO constraint to $600 N$ for this example, since the dynamics for optimal electrical power will be different." + "Next, we will create a new `PTO` object with this impedance matrix.\n", + "We will also update the definitions of our PTO constraint and additional dynamic forcing function to use the new object.\n", + "We will again set our maximum PTO force to 750 Newtons in this example, so we do not need to redefine `f_max`." ] }, { @@ -488,10 +692,9 @@ "pto_2 = wot.pto.PTO(ndof, kinematics, controller, pto_impedance_2, loss, name_2)\n", "\n", "## Update PTO constraints and forcing\n", - "f_max_2 = 600.0\n", "def const_f_pto_2(wec, x_wec, x_opt, waves):\n", " f = pto_2.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n", - " return f_max_2 - np.abs(f.flatten())\n", + " return f_max - np.abs(f.flatten())\n", "ineq_cons_2 = {'type': 'ineq', 'fun': const_f_pto_2}\n", "constraints_2 = [ineq_cons_2]\n", "f_add_2 = {'PTO': pto_2.force_on_wec}" @@ -501,7 +704,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Finally, let's update our `WEC` object with the new PTO constraint, then run our optimization problem. Note we're now using `average_power` instead of `mechanical_average_power` as our objective function." + "Finally, we will update our `WEC` object with the new PTO and objective function, then run our optimization problem.\n", + "Note we are now using `average_power` instead of `mechanical_average_power` as our objective function." ] }, { @@ -511,7 +715,6 @@ "outputs": [], "source": [ "# Update WEC\n", - "\n", "wec_2 = wot.WEC.from_bem(bem_data,\n", " constraints=constraints_2,\n", " friction=None,\n", @@ -538,15 +741,16 @@ "print(f'Optimal average electrical power: {opt_average_power} W')\n", "\n", "# Post-process\n", - "wec_fdom_2, wec_tdom_2 = wec_2.post_process(results_2[0], waves.sel(realization=0), nsubsteps)\n", - "pto_fdom_2, pto_tdom_2 = pto_2.post_process(wec_2, results_2[0], waves.sel(realization=0), nsubsteps)" + "wec_fdom_2, wec_tdom_2 = wec_2.post_process(wec_2, results_2, waves, nsubsteps)\n", + "pto_fdom_2, pto_tdom_2 = pto_2.post_process(wec_2, results_2, waves, nsubsteps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We will compare our optimal results to the unconstrained case to gain some insight into the effect of the constraint on the optimal PTO force. Let's do the same process as before, but unset the `constraints` parameter in a new `WEC` object." + "We will compare our optimal results to the unconstrained case to gain some insight into the effect of the constraint on the optimal PTO force.\n", + "Let's do the same process as before, but unset the `constraints` parameter in a new `WEC` object." ] }, { @@ -572,9 +776,9 @@ "opt_average_power = results_2_nocon[0].fun\n", "print(f'Optimal average electrical power: {opt_average_power} W')\n", "wec_fdom_2_nocon, wec_tdom_2_nocon = wec_2_nocon.post_process(\n", - " results_2_nocon[0], waves.sel(realization=0), nsubsteps)\n", + " wec_2_nocon, results_2_nocon, waves, nsubsteps)\n", "pto_fdom_2_nocon, pto_tdom_2_nocon = pto_2.post_process(\n", - " wec_2_nocon, results_2_nocon[0], waves.sel(realization=0), nsubsteps)" + " wec_2_nocon, results_2_nocon, waves, nsubsteps)" ] }, { @@ -582,7 +786,7 @@ "metadata": {}, "source": [ "Note that the optimal constrained PTO force follows the optimal unconstrained solution (sinusoidal) whenever the unconstrained solution is within the constraint. \n", - "When the constraint is active the optimal PTO force is the maximum PTO force of $600 N$." + "When the constraint is active the optimal PTO force is the maximum PTO force of 750 Newtons." ] }, { @@ -592,18 +796,21 @@ "outputs": [], "source": [ "plt.figure()\n", - "wec_tdom_2['pos'].plot(label='constrained')\n", - "wec_tdom_2_nocon['pos'].plot(label='unconstrained')\n", + "wec_tdom[0]['pos'].plot(label='optimal mechanical')\n", + "wec_tdom_2[0]['pos'].plot(label='optimal electrical, constrained')\n", + "wec_tdom_2_nocon[0]['pos'].plot(label='optimal electrical, unconstrained')\n", "plt.legend(loc='lower right')\n", "\n", "plt.figure()\n", - "pto_tdom_2['force'].plot(label='constrained')\n", - "pto_tdom_2_nocon['force'].plot(label='unconstrained')\n", + "pto_tdom[0]['force'].plot(label='optimal mechanical')\n", + "pto_tdom_2[0]['force'].plot(label='optimal electrical, constrained')\n", + "pto_tdom_2_nocon[0]['force'].plot(label='optimal electrical, unconstrained')\n", "plt.legend(loc='lower right')\n", "\n", "plt.figure()\n", - "pto_tdom_2['power'].loc['elec',:,:].plot(label='constrained')\n", - "pto_tdom_2_nocon['power'].loc['elec',:,:].plot(label='unconstrained')\n", + "pto_tdom[0]['power'].sel(type='mech').plot(label='optimal mechanical')\n", + "pto_tdom_2[0]['power'].sel(type='elec').plot(label='optimal electrical, constrained')\n", + "pto_tdom_2_nocon[0]['power'].sel(type='elec').plot(label='optimal electrical, unconstrained')\n", "plt.legend(loc='lower right')\n" ] }, @@ -611,7 +818,24 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The attentive user might have noticed that the amplitude of the position, force and power signals is about half the magnitude of the signals we plotted in the first part of the tutorial. We can see that optimizing for electrical power requires optimal state trajectories with smaller amplitudes. For most WECs the electrical power is the usable form of power, thus the WEC should be designed for electrical power and we can avoid over-designing, which would results from expecting the forces associated with the optimal trajectories for mechanical power maximisation." + "The attentive user might have noticed that the amplitude of the mechanical power is less compared to Part 1 of the tutorial.\n", + "We can see that optimizing for electrical power requires optimal state trajectories with less reactive mechanical power (i.e. power that is put into the system).\n", + "\n", + "The PTO force trajectory for optimizing mechanical power is saturated at the maximum for longer compared to the electrical power.\n", + "This could inform the WEC designer optimizing for mechanical power to consider larger components that would not be utilized at their limit as frequently.\n", + "However, the electrical power (_not_ the mechanical power) is the usable form of power, thus designing the WEC for optimal electrical power does not indicate a need for larger components and prevents this over-design.\n", + "\n", + "The Sankey power flow diagram confirm this observation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p_flows_2 = wot.utilities.calculate_power_flows(wec_2, pto_2, results_2, waves, intrinsic_impedance)\n", + "wot.utilities.plot_power_flow(p_flows_2)" ] }, { @@ -619,16 +843,18 @@ "metadata": {}, "source": [ "## 3. Control co-design of the WEC geometry for maximum electrical power\n", - "The first two examples only used the inner optimization loop in WecOptTool to optimize PTO power. Here in Part 3 we bring it all together and show how to use both the inner and outer optimization loops in WecOptTool to do control co-optimization of a hull design in conjunction with an optimal controller for electrical power.\n", + "The first two examples only used the inner optimization loop in WecOptTool to optimize PTO power.\n", + "In Part 3, we bring it all together and show how to use both the inner and outer optimization loops in WecOptTool to do control co-optimization of a hull design in conjunction with an optimal controller for electrical power.\n", + "\n", "Again, we use the WaveBot WEC in one degree of freedom in regular waves. \n", "The goal is to **find the optimal keel radius** (`r2`) that maximizes the average produced electrical power, while maintaining a constant hull volume. \n", - "A constant volume is achieved by setting the height of the conical section (`h2`) in conjunction with the keel radius (`r2`).\n", + "A constant volume is achieved by setting the height of the conical section (`h2`) in conjunction with the keel radius.\n", "\n", - "This example demonstrates a complete case of the types of optimization studies WecOptTool is meant for. \n", - "The main optimization (outer optimization loop) is to find the optimal geometry (radius `r2`), and for each geometry considered the optimal PTO force (inner optimization loop) will be found.\n", - "The inner loop was showcased in Example 2 and uses a gradient-based optimization method, with the gradients obtained with automatic differentiation. \n", - "The outer loop optimization is for the user to setup. \n", - "In this example, we will do a simple *brute force* optimization using `scipy.optimize.brute`. \n", + "This example demonstrates a complete case of the control co-optimization studies WecOptTool is meant for. \n", + "The inner optimization loop finds the control trajectory that produces the optimal PTO force for a given hull geometry, and the outer optimization loop finds the optimal hull geometry _amongst designs with optimal control trajectories_.\n", + "\n", + "The inner loop is consolidated into the `WEC.solve()` method, but the outer loop needs to be configured by the user for their particular design problem.\n", + "In this example, we will do a simple *brute force* optimization using `scipy.optimize.brute` (click [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brute.html) for documentation). \n", "\n", "![Device Diagram](https://live.staticflickr.com/65535/51751577441_515afec334_z.jpg) \n", "
\n", @@ -681,7 +907,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next we will define an objective function for our design optimization problem. We use the same workflow illustrated in Part 2 to set up a WaveBot device and solve for the optimal solution, but wrap this in a function definition which can set `r2` and (indirectly) `h2`." + "Next we will define an objective function for our design optimization problem.\n", + "We use the same workflow used in Part 2 to set up a WaveBot device and solve for the optimal solution, but wrap this in a function definition which sets `r2` and (indirectly) `h2`." ] }, { @@ -744,7 +971,7 @@ "\n", " # Set PTO constraint and additional dynamic force\n", " nsubsteps = 4\n", - " f_max = 600.0\n", + " f_max = 750.0\n", "\n", " def const_f_pto(wec, x_wec, x_opt, waves):\n", " f = pto.force_on_wec(wec, x_wec, x_opt, waves, nsubsteps)\n", @@ -762,12 +989,6 @@ " f_add=f_add,\n", " )\n", "\n", - " # Waves\n", - " wfreq = 0.3\n", - " amplitude = 0.0625\n", - " phase = -40\n", - " waves_3 = wot.waves.regular_wave(f1, nfreq, wfreq, amplitude, phase)\n", - "\n", " # Objective function\n", " obj_fun = pto.average_power\n", " nstate_opt = 2*nfreq\n", @@ -792,12 +1013,11 @@ "metadata": {}, "source": [ "### Solve\n", - "Finally, we may call this objective function with an optimization algorithm. \n", - "Here, a simple *brute force* optimization approach is used for illustrative purposes, but any variety of options could be applied. \n", + "Finally, we call our wrapped function using `scipy.optimize.brute`.\n", "The optimization algorithm will call our objective function, which in turn will create a new WaveBot hull, run the necessary BEM calculations for the hull, and find the PTO force that provides the most electric power for that hull. \n", "This process will be conducted for the range of `r2` values that we specify.\n", "\n", - "_(note: the cell below will likely take 5+ minutes to run on a standard personal computer)_" + "_(note: the cell below will take longer to run than the cells above, up to several minutes)_" ] }, { @@ -820,8 +1040,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`." ] }, { @@ -845,7 +1064,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that in this case the magnitude of average power between the different keel radii is rather small, this is because the PTO force constraint is active most of the time, therefore all considered geometries perform similarily. If you remove the PTO constraint and re-run the co-optimization study you will see that the impact of radius on average electrical power is significantly higher." + "Note that in this case the difference in average power between the different keel radii is rather small.\n", + "This is because the PTO force constraint is active most of the time, therefore all considered geometries perform similarily.\n", + "If you remove the PTO constraint and re-run the co-optimization study, you will see that the impact of radius on average electrical power is significantly higher." ] } ], @@ -865,7 +1086,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.11.4" }, "vscode": { "interpreter": { diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb index 83df7e2d3..8d42ede4e 100644 --- a/examples/tutorial_2_AquaHarmonics.ipynb +++ b/examples/tutorial_2_AquaHarmonics.ipynb @@ -5,13 +5,13 @@ "metadata": {}, "source": [ "# Tutorial 2 - AquaHarmonics\n", - "The goal of this tutorial is to illustrate a more realistic model of a PTO, including non-linear power conversion chain. \n", + "The goal of this tutorial is to illustrate a more realistic model of a PTO, including nonlinear power conversion chain. \n", "It uses the [AquaHarmonics](https://aquaharmonics.com/wec_vis/) device in one degree of freedom in regular waves, models the PTO generator using a motor power loss map and adds realistic constraints, including generator maximum torque and min/max line tension.\n", "\n", "This tutorial is comprises two parts, where the second section builds upon the first.\n", "\n", "1. [Optimal control with a loss map](#1.-Optimal-control-with-a-loss-map)\n", - "2. [Control co-design of line pre-tension/ballast with non-linear power conversion](#2.-Control-co-design-of-line-pre-tension/ballast)\n", + "2. [Control co-design of line pretension/ballast with non-linear power conversion](#2.-Control-co-design-of-line-pretension/ballast)\n", "\n", "![AquaHarmonics device](https://aquaharmonics.com/wec_vis.png)" ] @@ -27,11 +27,12 @@ "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" + "import wecopttool as wot\n", + "\n", + "## set colorblind-friendly colormap for plots\n", + "plt.style.use('tableau-colorblind10')\n", + "cc = plt.rcParams['axes.prop_cycle'].by_key()['color']" ] }, { @@ -41,12 +42,13 @@ "## 1. Optimal control with a loss map\n", "\n", "### WEC object\n", - "We will go through a number of steps to create a `WEC` object, with which we can solve for the optimal control trajectory.\n", - "To compose the `WEC` object, we will need to define the mesh, degrees of freedom, mass and hydrostatic properties, linear hydrodynamic coefficients (from a BEM solution), any additional dynamic forces (e.g. PTO force), and constraints (e.g. minimum line tension). \n", + "We will go through a number of steps to create a `WEC` object, with we can use to solve for the optimal control trajectory.\n", + "See Tutorial 1 for detailed instructions on what is needed to create the `WEC` object.\n", "\n", "#### Geometry\n", "First we create a surface mesh for the WEC hull, and quickly visualize the cross-sectional shape of the hull.\n", - "The AquaHarmonics hull mesh is already included with WecOptTool—take a look at the `geometry.py` module if you're interested in seeing how it was created.\n", + "The AquaHarmonics hull mesh is already included with WecOptTool, much like the WaveBot mesh in Tutorial 1.\n", + "Take a look at the `geom.py` module [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) if you are interested in seeing how it was created.\n", "Note that the lower cylindrical section of the hull is open to the sea." ] }, @@ -87,7 +89,8 @@ "source": [ "#### Hydrostatics and mass\n", "The AquaHarmonics device is positively buoyant (i.e., the buoyancy force is greater than the force due to gravity).\n", - "Therefore, we'll set the rigid-body mass manually, but allow the hydrostatic stiffness to be set automatically by run_bem() based on the mesh. We will also calculate the displaced volume and mass from the mesh (before manually defining the mass of the FloatingBody), which we will need later for defining forces and constraints." + "Therefore we will set the rigid-body mass manually (unlike Tutorial 1), but allow the hydrostatic stiffness to be set automatically by the `run_bem` function based on the mesh.\n", + "We will also calculate the displaced volume and mass from the mesh (before manually defining the mass of the FloatingBody), which we will need later for defining forces and constraints." ] }, { @@ -100,10 +103,36 @@ "rho = 1025\n", "fb.center_of_mass = [0, 0, 0]\n", "fb.rotation_center = fb.center_of_mass\n", - "displaced_mass = fb.compute_rigid_body_inertia(rho=rho).values # [kg]\n", - "displacement = displaced_mass/rho # [m^3] \n", + "displaced_mass = fb.compute_rigid_body_inertia(rho=rho).values # kg\n", + "displacement = displaced_mass/rho # m^3\n", "\n", - "fb.mass = np.atleast_2d(5e3) # [kg]" + "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, `f1`, and number of frequencies, `nfreq`, to ensure that 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)" ] }, { @@ -113,8 +142,8 @@ "#### 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", - "For `WecOptTool`, these must be a regularly spaced frequency array." + "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 +152,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)" @@ -143,23 +170,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axes = plt.subplots(3,1)\n", - "bem_data['added_mass'].plot(ax = axes[0])\n", - "bem_data['radiation_damping'].plot(ax = axes[1])\n", - "axes[2].plot(bem_data['omega'],np.abs(np.squeeze(bem_data['diffraction_force'].values)), color = 'orange')\n", - "axes[2].set_ylabel('abs(diffraction_force)', color = 'orange')\n", - "axes[2].tick_params(axis ='y', labelcolor = 'orange')\n", - "ax2r = axes[2].twinx()\n", - "ax2r.plot(bem_data['omega'], np.abs(np.squeeze(bem_data['Froude_Krylov_force'].values)), color = 'blue')\n", - "ax2r.set_ylabel('abs(FK_force)', color = 'blue')\n", - "ax2r.tick_params(axis ='y', labelcolor = 'blue')\n", - "\n", - "for axi in axes:\n", - " axi.set_title('')\n", - " axi.label_outer()\n", - " axi.grid()\n", - " \n", - "axes[-1].set_xlabel('Frequency [rad/s]')" + "[(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)] = wot.utilities.plot_hydrodynamic_coefficients(bem_data)" ] }, { @@ -233,7 +244,7 @@ "source": [ "##### Generator loss map\n", "The generator used by the AquaHarmonics devices is a motor from a Nissan Leaf.\n", - "`WecOptTool` currently uses loss maps, as efficiency maps are unable to account for losses that occur when the mechanical power is zero (e.g., when the torque is nonzero, but the shaft speed is zero).\n", + "WecOptTool currently uses loss maps, as efficiency maps are unable to account for losses that occur when the mechanical power is zero (e.g., when the torque is nonzero, but the shaft speed is zero).\n", "Thus, we use an [approximate model for the losses](https://en.wikipedia.org/wiki/Joule_heating) in this example." ] }, @@ -268,9 +279,9 @@ "metadata": {}, "outputs": [], "source": [ - "rot_max = 10000*2*np.pi/60 # [rad/s]\n", - "torque_max = 300 # [Nm]\n", - "power_max = 80e3 # [W]" + "rot_max = 10000*2*np.pi/60 # rad/s\n", + "torque_max = 300 # Nm\n", + "power_max = 80e3 # W" ] }, { @@ -278,7 +289,7 @@ "metadata": {}, "source": [ "Finally, we can plot the motor power loss surface to see how it looks.\n", - "Here, the independent variables are the motor speed in [rad/s] and the motor torque in [Nm].\n", + "Here, the independent variables are the motor speed (in rad/s) and the motor torque (in Nm).\n", "Note that we limit the plot by the power limit." ] }, @@ -298,13 +309,12 @@ "ax = [fig.add_subplot(1, 2, 1, projection=\"3d\"),\n", " fig.add_subplot(1, 2, 2)]\n", "\n", - "ax[0].plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0)\n", + "ax[0].plot_surface(X, Y, Z, cmap=cm.viridis, linewidth=0)\n", "ax[0].set_xlabel('Shaft speed [rad/s]')\n", "ax[0].set_ylabel('Torque [Nm]')\n", "ax[0].set_zlabel('Power loss [kW]')\n", - "# ax[0].set_zlim([0, 1.1])\n", "\n", - "contour = ax[1].contourf(X, Y, Z, cmap=cm.coolwarm)\n", + "contour = ax[1].contourf(X, Y, Z, cmap=cm.viridis)\n", "plt.colorbar(contour, label=\"Power loss [kW]\")\n", "ax[1].set_xlabel('Shaft speed [rad/s]')\n", "ax[1].set_ylabel('Torque [Nm]')\n", @@ -344,10 +354,10 @@ "metadata": {}, "source": [ "#### Additional Forces\n", - "In WecOptTool, hydrodynamic forces are considered automatically.\n", - "Thus to account for other phenomena we create ''**additional forces**'':\n", + "Beyond hydrodynamic loading, this model needs to account for these additional phenomena:\n", "\n", - "* **Buoyancy, gravity and pretension** - While the buoyancy and gravity phenomena are often lumped together, we will keep them separate here. This is useful because the AquaHarmonics devices has a pretension force, so the buoyancy and gravity forces are not balanced at the neutral position and there is a pretension force from the tether/air spring.\n", + "* **Buoyancy, gravity and pretension** - While buoyancy and gravity are often lumped together in similar models, we will keep them separate here.\n", + "This is useful because the AquaHarmonics devices has a pretension force, so the buoyancy and gravity forces are not balanced at the neutral position and there is a pretension force from the tether/air spring.\n", "* **Passive PTO forces** - Inertia forces due to finite rigid body inertia of gears and friction forces within the drive train are captured by this function.\n", "* **PTO line force** - The total force imposed on the WEC from its tether is the sum of the PTO force due to action by the generator and the pretension force from the air spring.\n", "The generator torque will be our control state for which we will derive the optimal trajectory." @@ -408,13 +418,16 @@ "Each of these will be implemented as inequality constraints that will be passed to the numerical optimization solver.\n", "\n", "Note that we are imposing the constraints at more points than the dynamics by setting `nsubsteps>1`.\n", - "This ensures that the constraint is properly maintained (see the Theory page in the WecOptTool documentation for more details).\n", + "This ensures that the constraint is properly maintained (see the [Theory](https://sandialabs.github.io/WecOptTool/theory.html) documentation for more details).\n", "\n", " * **Peak torque** - The absolute motor torque cannot exceed this value.\n", " * **Maximum rotational speed** - The absolute motor speed cannot exceed this value.\n", " * **Maximum power** - The maximum mechanical power (i.e., product of torque and velocity) cannot exceed this value.\n", - " * **Minumum line tension** - In addition to these limitations on the hardware, we will also constrain the solution to prevent the tether tension from going below a certain threshold. We absolutely cannot allow the line tension to be less than zero, or else it will go slack. In practice, it is prudent to set some nonzero limit for the tether tension.\n", - " * **Mean torque** - While solutions with an nonzero mean positions may indeed be valid, we want to rely on a single calculation of the linear hydrodynamics at a mean position of 0. To avoid violating this approach, we will constrain the mean WEC position to zero." + " * **Minimum line tension** - In addition to these limitations on the hardware, we will also constrain the solution to prevent the tether tension from going below a certain threshold.\n", + " We absolutely cannot allow the line tension to be less than zero, or else it will go slack.\n", + " In practice, it is prudent to set some nonzero limit for the tether tension.\n", + " * **Mean torque** - While solutions with an nonzero mean positions may indeed be valid, we want to rely on a single calculation of the linear hydrodynamics at a mean position of 0.\n", + " To avoid violating this approach, we will constrain the mean WEC position to zero." ] }, { @@ -424,9 +437,10 @@ "outputs": [], "source": [ "# Generator constraints\n", - "torque_peak_max = 280 #[Nm] \n", - "rot_speed_max = 10000*2*np.pi/60 #[rad/s]\n", - "# mooring/PTO line constraint\n", + "torque_peak_max = 280 # Nm\n", + "rot_speed_max = 10000*2*np.pi/60 # rad/s\n", + "\n", + "# Mooring/PTO line constraint\n", "min_line_tension = -1000\n", "\n", "nsubsteps = 2\n", @@ -467,9 +481,9 @@ "metadata": {}, "source": [ "#### WEC object\n", - "Finally, we can use all the different components we've developed thus far to construct a `WEC` object:\n", + "Finally, we can use all the different components developed thus far to construct the `WEC` object:\n", "\n", - " * **BEM data** - defines the hydrodynamics of the hull and includes hydrostatics\n", + " * **BEM data** - defines the hydrostatics and hydrodynamics of the hull\n", " * **Constraints** - limitations on the hardware (max power, max torque, etc.) and the constraint to prevent the tether from going slack\n", " * **Additional forces** - this captures all of the forces on the WEC that are not due to the interaction between the hull and water (PTO, etc.)" ] @@ -489,28 +503,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": {}, @@ -573,8 +565,7 @@ "metadata": {}, "source": [ "### Post-process and plotting\n", - "Next, we post-process the results to allow us to more easily visualize them in a series of plots.\n", - "Note that because the `WEC` and `PTO` objects are distinct (the WEC object really only knows what the force from the PTO is, not how that force is obtained), we create two separate results objects (in each case, we get results in the frequency domain and time domain)." + "Next, we post-process the results to allow us to more easily visualize them in a series of plots." ] }, { @@ -583,15 +574,15 @@ "metadata": {}, "outputs": [], "source": [ - "wec_fdom, wec_tdom = wec.post_process(results[0], waves.sel(realization=0), nsubsteps=nsubsteps)\n", - "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves.sel(realization=0), nsubsteps=nsubsteps)" + "wec_fdom, wec_tdom = wec.post_process(wec, results, waves, nsubsteps=nsubsteps)\n", + "pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can inspect each of these [Xarray Datasets](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html) for more details.\n", + "We can inspect each list of `xarray.Dataset` objects for more details.\n", "As an example the post-processed WEC time-domain results include the following:" ] }, @@ -627,10 +618,10 @@ "\n", "# Excitation and velocity\n", "ax1 = ax[0].twinx()\n", - "force_excitation = wec_tdom.force.sel(type=['Froude_Krylov', 'diffraction'])\n", + "force_excitation = wec_tdom[0].force.sel(type=['Froude_Krylov', 'diffraction'])\n", "force_excitation = force_excitation.sum('type')/1e3\n", "force_excitation.plot(ax=ax[0], color='k')\n", - "wec_tdom.vel.plot(ax=ax1, color='orange')\n", + "wec_tdom[0].vel.plot(ax=ax1, color='orange')\n", "ax1.set_ylabel('Velocity [m/s]', color='orange')\n", "ax1.tick_params(axis='y', labelcolor='orange')\n", "ax1.set_title('')\n", @@ -638,7 +629,7 @@ "ax[0].set_ylabel('Excitation force [kN]')\n", "\n", "# Forces\n", - "(wec_tdom.force.sel(type='PTO')/1e3).plot(\n", + "(wec_tdom[0].force.sel(type='PTO')/1e3).plot(\n", " ax=ax[1], \n", " label='PTO force in WEC frame',\n", ")\n", @@ -661,21 +652,21 @@ "ax[1].legend(loc=1)\n", "\n", "# Torque\n", - "pto_tdom.force.plot(\n", + "pto_tdom[0].force.plot(\n", " ax=ax[2], \n", " linestyle='solid',\n", " label='PTO torque in PTO frame',\n", ")\n", "ax[2].plot(\n", - " pto_tdom.time, \n", - " 1*torque_peak_max * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, \n", + " 1*torque_peak_max * np.ones(pto_tdom[0].time.shape),\n", " color='black', \n", " linestyle='dotted',\n", " label=f'Peak torque limit ($\\pm${torque_peak_max}Nm)',\n", ")\n", "ax[2].plot(\n", - " pto_tdom.time, \n", - " -1*torque_peak_max * np.ones(pto_tdom.time.shape), \n", + " pto_tdom[0].time, \n", + " -1*torque_peak_max * np.ones(pto_tdom[0].time.shape), \n", " color='black', \n", " linestyle='dotted',\n", ")\n", @@ -683,12 +674,9 @@ "ax[2].legend(loc='upper right',)\n", "\n", "# Power\n", - "(pto_tdom['power'].loc['mech',:,:]/1e3).plot(ax=ax[3], label='Mech. power')\n", - "(pto_tdom['power'].loc['elec',:,:]/1e3).plot(ax=ax[3], linestyle='dashed', \n", + "(pto_tdom[0]['power'].sel(type='mech')/1e3).plot(ax=ax[3], label='Mech. power')\n", + "(pto_tdom[0]['power'].sel(type='elec')/1e3).plot(ax=ax[3], linestyle='dashed', \n", " label=\"Elec. power\")\n", - "# ax[3].axhline(-1*power_max/1e3, linestyle='dotted', color='black', \n", - " # label=f'Max mech. power ($\\pm${power_max/1e3:.0f}kW)')\n", - "# ax[3].axhline(1*power_max/1e3, linestyle='dotted', color='black')\n", "ax[3].grid(color='0.75', linestyle='-', linewidth=0.5, axis='x')\n", "ax[3].legend(loc='upper right')\n", "ax[3].set_title('')\n", @@ -698,7 +686,7 @@ " axi.set_title('')\n", " axi.grid()\n", " axi.label_outer()\n", - " axi.autoscale(axis='x', tight=True)\n" + " axi.autoscale(axis='x', tight=True)" ] }, { @@ -717,35 +705,37 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(8,6))\n", - "ax.plot(pto_tdom.vel, pto_tdom.force, color='k')\n", - "contour = ax.contourf(X, Y, Z, levels=20, cmap=cm.coolwarm)\n", + "ax.plot(pto_tdom[0].vel, pto_tdom[0].force, color=cc[1], label='Optimal control trajectory')\n", + "contour = ax.contourf(X, Y, Z, levels=20, cmap=cm.viridis)\n", "plt.colorbar(contour, label=\"Power loss [kW]\")\n", "ax.grid(which='major', linestyle='--')\n", "ax.axvline(0, color='k', lw=1)\n", "ax.axhline(0, color='k', lw=1)\n", "ax.set_ylim([-torque_max, torque_max])\n", "ax.set_ylabel('Torque [Nm]')\n", - "ax.set_xlabel('Shaft speed [rad/s]')" + "ax.set_xlabel('Shaft speed [rad/s]')\n", + "fig.legend(loc='lower right')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Control co-design of line pre-tension/ballast\n", + "## 2. Control co-design of line pretension/ballast\n", "[Part 1](#1.-Optimal-control-with-a-loss-map) of this tutorial set up the model and solved for the optimal control trajectory of a single fixed device design.\n", "We now will do control co-design to answer a real, practical question the designers have about the AquaHarmonics design:\n", - "**How big of an effect does the mass vs. line pre-tension have on the output power?** \n", + "**How big of an effect does the mass vs. line pretension have on the output power?** \n", "\n", - "To do this study we will define the maximum mass as the mass that results in the pre-tension being equal to the minimum pretension, and define that as a mass ratio of `1`.\n", + "To do this study we will define the maximum mass as the mass that results in the pretension being equal to the minimum pretension, and define that as a mass ratio of 1.\n", "Note that this maximum mass is slightly smaller than the displaced mass, in order to maintain some positive net buoyancy and thus achieve the minimum line tension.\n", "We will then consider different designs consisting of different mass ratios and see how the output power varies. \n", - "For each design the optimal controller for that design will be found (as in [Part 1](#1.-Optimal-control-with-a-loss-map))\n", + "For each design the optimal controller for that design will be found, as in Part 1.\n", "\n", "### Problem setup\n", - "The first step is to create a function that encapsulates solving for the optimal control trajectory to maximize power (i.e., [Part 1](#1.-Optimal-control-with-a-loss-map) of this tutorial) as an inner loop nested within an outer loop that varies the mass ratio.\n", + "The first step is to create a function that encapsulates solving for the optimal control trajectory to maximize power (i.e., Part 1 of this tutorial) to nest within an outer optimization loop.\n", + "This is similar to part 3 of Tutorial 1, except we now vary the mass ratio instead of the WEC hull radius.\n", "To accomplish this, we will create a function which takes the mass fraction as an input and returns the average electrical power as an output.\n", - "We consider the same regular wave operational condition as in [Part 1](#1.-Optimal-control-with-a-loss-map). " + "We consider the same regular wave operational condition used in Part 1." ] }, { @@ -841,9 +831,9 @@ "In this wave, the most power is generated when the mass fraction is 0.7.\n", "This is a scenario when the ballast is such that the total rigid body mass is 70% of the maximum, which allows for a pretension that is approximately 21 times the minimum pretension.\n", "\n", - "As the mass fraction is increased, pre-tension approaches the minimum line tension.\n", + "As the mass fraction is increased, pretension approaches the minimum line tension.\n", "At the limit (mass fraction of unity), the WEC can produce no power, as any action by the motor would slack the mooring line (recall that we forced the solution be symmetric).\n", - "Recall from the results in [Part 1](#1.-Optimal-control-with-a-loss-map) that we observed the motor torque saturating to avoid violating the minimum line tension constraint.\n", + "Recall from [the results in Part 1](#time-histories) that we observed the motor torque saturating to avoid violating the minimum line tension constraint.\n", "The red curve corresponding to the y-axis on the right-hand side of the plot shows the minimum absolute line tension.\n", "We can see that for mass fractions equal to and greater than 0.7 the line tension reaches the minimum value, thus bringing that constraint into effect and limiting the generated power.\n", "\n", @@ -863,11 +853,13 @@ "metadata": {}, "outputs": [], "source": [ + "cc = plt.rcParams['axes.prop_cycle'].by_key()['color']\n", + "\n", "fig, ax = plt.subplots(figsize=(8,6))\n", - "ax.plot(res[2], res[3]/1e3, 'bo-')\n", + "ax.plot(res[2], res[3]/1e3, 'o-', color=cc[0])\n", "ax.set_xlabel(\"Mass ratio, $m/m_{max}$\")\n", - "ax.set_ylabel(\"Average electric power, $\\overline{P}_e$ [kW]\", color=\"b\")\n", - "ax.tick_params(axis='y', labelcolor=\"b\")\n", + "ax.set_ylabel(\"Average electric power, $\\overline{P}_e$ [kW]\", color=cc[0])\n", + "ax.tick_params(axis='y', labelcolor=cc[0])\n", "\n", "# second tick\n", "new_tick_locations = [0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]\n", @@ -879,24 +871,17 @@ "ax2.set_xlim(ax.get_xlim())\n", "ax2.set_xticks(new_tick_locations)\n", "ax2.set_xticklabels(tick_function(new_tick_locations))\n", - "ax2.set_xlabel(\"Pre-tension ratio, $t/t_{min}$\")\n", + "ax2.set_xlabel(\"Pretension ratio, $t/t_{min}$\")\n", "\n", "ax3 = ax.twinx()\n", "ax3.axhline(-1*min_line_tension/1e3, linestyle='dotted', color='black',\n", " label='Min. line tension')\n", - "ax3.plot(res[2], np.abs(max_torque)/1e3, 'ro-')\n", - "ax3.set_ylabel(\"Min. line tension [kN]\", color=\"r\")\n", - "ax3.tick_params(axis='y', labelcolor=\"r\")\n", + "ax3.plot(res[2], np.abs(max_torque)/1e3, 'o-', color=cc[1])\n", + "ax3.set_ylabel(\"Min. line tension [kN]\", color=cc[1])\n", + "ax3.tick_params(axis='y', labelcolor=cc[1])\n", "ax3.set_yticks([1, 5, 10, 15, 20, 25])\n", "ax3.legend(loc=9)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -915,7 +900,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.11.4" }, "vscode": { "interpreter": { diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb index 065ae977f..7aee37448 100644 --- a/examples/tutorial_3_LUPA.ipynb +++ b/examples/tutorial_3_LUPA.ipynb @@ -6,16 +6,20 @@ "source": [ "# Tutorial 3 - LUPA\n", "\n", - "The main goal of this tutorial is to demonstrate how to set up a WEC design problem with a more complex and realistic setup. We use the [Lab Upgrade Point Absorber (LUPA)](https://pmec-osu.github.io/LUPA/) device, an open source two-body heaving point absorber under development by Oregon State University. A deep dive video demonstration of the LUPA device and its features can be viewed [here](https://www.youtube.com/watch?v=gCcAu7H9lQI). We will numerically replicate LUPA testing in the Large Wave Flume at the [O.H. Hinsdale Wave Research Laboratory](https://engineering.oregonstate.edu/wave-lab) in order to provide further design optimization for the WEC device concept. This tutorial builds on the previous ones and introduces:\n", + "The main goal of this tutorial is to demonstrate how to set up a WEC design problem with a more complex and realistic setup.\n", + "We use the [Lab Upgrade Point Absorber (LUPA)](https://pmec-osu.github.io/LUPA/) device, an open source two-body heaving point absorber under development by Oregon State University.\n", + "A deep dive video demonstration of the LUPA device and its features can be viewed [here](https://www.youtube.com/watch?v=gCcAu7H9lQI).\n", + "We will numerically replicate LUPA testing in the Large Wave Flume at the [O.H. Hinsdale Wave Research Laboratory](https://engineering.oregonstate.edu/wave-lab) in order to provide further design optimization for the WEC device concept.\n", + "This tutorial builds on the previous ones by introducing:\n", " \n", - "- a WEC comprised of multiple bodies\n", - "- setting up a WEC with multiple degrees of freedom (DOF) using generalized modes\n", - "- more complex PTO kinematics that depend on more than one WEC DOF \n", - "- realistic constraints including generator maximum and continuous torque, and maximum rotational speed\n", - "- irregular waves\n", - "- mooring system dynamics\n", + "* a WEC comprised of multiple bodies\n", + "* setting up a WEC with multiple degrees of freedom (DOF) using generalized modes\n", + "* more complex PTO kinematics that depend on more than one WEC DOF \n", + "* realistic constraints including generator maximum and continuous torque, and maximum rotational speed\n", + "* irregular waves\n", + "* mooring system dynamics\n", "\n", - "A secondary goal is for this notebook to serve as a tool for those who are planning to run experiments with the LUPA device, to inform their experiment design. \n", + "A secondary goal for this tutorial is to serve as a tool for those who are planning to run experiments with the LUPA device, to inform their experiment design. \n", "\n", "As with previous tutorials, this tutorial consists of two parts, with the second section building upon the first.\n", "\n", @@ -32,6 +36,7 @@ "outputs": [], "source": [ "import os\n", + "\n", "import gmsh, pygmsh\n", "import capytaine as cpy\n", "import autograd.numpy as np\n", @@ -52,9 +57,18 @@ "## 1. Optimal control of a two-body WEC\n", "\n", "### WEC geometry\n", - "The creation of the `WEC` object is fundamentally identical to previous tutorials, where we use meshes of the WEC to create Capytaine `FloatingBody` objects, run BEM using Capytaine, and use the `WEC.from_bem()` method to create the `WEC` object. The key here is that the LUPA is a two-body device (consisting of a float and a spar), which move independently in heave but in unison for all other degrees of freedom. To model this in WecOptTool, we can create a `FloatingBody` object for each body separately with a heave DOF, combine them into a single object afterwards, and be sure the combined mass and inertia properties are properly set.\n", + "The creation of the `WEC` object is fundamentally identical to previous tutorials, where we use meshes of the WEC to create Capytaine `FloatingBody` objects, run BEM using Capytaine, and use the `WEC.from_bem()` method to create the `WEC` object.\n", + "The key here is that the LUPA is a two-body device (consisting of a float and a spar), which move independently in heave but in unison for all other degrees of freedom.\n", + "To model this in WecOptTool, we can create a `FloatingBody` object for each body separately with a heave DOF, combine them into a single object afterwards, and be sure the combined mass and inertia properties are properly set.\n", + "\n", + "We will analyze the device in its four planar degrees of freedom:\n", + "\n", + "* Heave of the buoy\n", + "* Heave of the spar\n", + "* Combined device surge\n", + "* Combined device pitch\n", "\n", - "We will analyze the device in its four planar degrees of freedom: the heave of the buoy, heave of the spar, combined device surge, and combined device pitch, for a total of 4 degrees of freedom. Here we are using the generalized modes approach. Alternatively we could include all 3 planar DOF for each body separately and add two constraints for the pitch and surge to be equal for both bodies. \n", + "Here we are using the generalized modes approach; an alternative solution would be to include all 3 planar DOF for each body separately and add two constraints for the pitch and surge to be equal for both bodies.\n", "\n", "#### LUPA properties\n", "The mass properties of the LUPA have been provided from measurements of the physical device by Oregon State University, as follows:" @@ -86,7 +100,12 @@ "source": [ "#### Mesh creation of the float\n", "\n", - "Here we create the mesh based on the dimensions provided by Oregon State University using *pygmsh* as in the other tutorials. The float has a hole where the spar passes through it; OSU has found that this hole creates a large spike in the BEM results at about 4.5 rad/s. This is not much energy in the waves at this frequency in the wave spectrum we will be using, so this spike should not significantly affect our results. If you would like to remove this spike and create smoother BEM results, set `hole = False` in the cell below." + "Here we create the mesh based on the dimensions provided by Oregon State University using `pygmsh`, available [here](https://pypi.org/project/pygmsh/).\n", + "This is the same package used by `geom.py` (click [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) for API documentation) containing the predefined WaveBot and AquaHarmonics meshes used in the previous tutorials.\n", + "\n", + "The float has a hole where the spar passes through it, though OSU has found that this hole creates a large spike in the BEM results at about 4.5 rad/s.\n", + "There is not much energy in the waves at this frequency in the wave spectrum we will be using, so this spike should not significantly affect our results.\n", + "If you would like to remove this spike and create smoother BEM results, set `hole = False` in the cell below." ] }, { @@ -102,7 +121,7 @@ "h1 = 0.5 \n", "h2 = 0.21\n", "freeboard = 0.3\n", - "r3 = 0.10/2 + 0.05 # hole radius # TODO\n", + "r3 = 0.10/2 + 0.05 # hole radius\n", "hole = True # set to False to remove spike in BEM results\n", "with pygmsh.occ.Geometry() as geom:\n", " gmsh.option.setNumber('Mesh.MeshSizeFactor', mesh_size_factor)\n", @@ -123,7 +142,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Again, we will only add the heave DOF for now. The surge and pitch will be added after we combine the two `FloatingBody` objects." + "Again, we will only add the heave DOF for now.\n", + "The surge and pitch will be added after we combine the two `FloatingBody` objects." ] }, { @@ -141,7 +161,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can now visualize the mesh. " + "We can now visualize the mesh:" ] }, { @@ -232,7 +252,11 @@ "source": [ "#### Combined `FloatingBody`\n", "\n", - "With both WEC bodies defined separately, we should first define the respective centers of mass and rotation centers for the `FloatingBody`s. Then, we can create a union of the bodies and define the properties for the overall LUPA device. At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension. The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia). This is also when we can specify the surge and pitch degrees of freedom.\n", + "With both WEC bodies defined separately, we can now define the respective centers of mass and rotation centers of the bodies.\n", + "We will then create a union of the bodies and define the properties for the overall LUPA device.\n", + "At the equilibrium position the float is neutrally buoyant while the spar is positively buoyant and requires mooring pre-tension.\n", + "The combined center of mass and moment of inertia can be found by using the given values weighted by the mass (via the parallel axis theorem for the moment of inertia).\n", + "This is also when we can specify the surge and pitch degrees of freedom.\n", "\n", "We are using the density of *fresh* water, $\\rho = 1000 kg/m^3$, since we are modeling LUPA in a wave flume." ] @@ -247,16 +271,16 @@ "rho = 1000\n", "\n", "# mass properties float\n", - "mass_float = float_mass_properties['mass'] \n", + "mass_float = float_mass_properties['mass']\n", "cm_float = np.array(float_mass_properties['CG'])\n", - "pitch_inertia_float = float_mass_properties['MOI'][1] \n", + "pitch_inertia_float = float_mass_properties['MOI'][1]\n", "float_fb.center_of_mass = cm_float\n", "float_fb.rotation_center = float_fb.center_of_mass\n", "\n", "# mass properties spar\n", - "mass_spar = spar_mass_properties['mass'] \n", + "mass_spar = spar_mass_properties['mass']\n", "cm_spar = np.array(spar_mass_properties['CG'])\n", - "pitch_inertia_spar = spar_mass_properties['MOI'][1] \n", + "pitch_inertia_spar = spar_mass_properties['MOI'][1]\n", "spar_fb.center_of_mass = cm_spar\n", "spar_fb.rotation_center = spar_fb.center_of_mass\n", "\n", @@ -288,8 +312,11 @@ "metadata": {}, "source": [ "#### Define Hydrostatics Manually\n", - "\n", - "When combining multiple bodies, we need to be careful to set up the hydrostatics correctly. The bodies move separately in heave but move together in surge and pitch. Therefore, we should define the individual heave inertia values for each body and define the total surge and pitch inertia values. The inertia then needs to be reformatted as a DataArray for Capytaine. Lastly, the hydrostatic stiffness can be calculated for the total immersed body." + "We need to be carefully set up the hydrostatic correctly when combining multiple bodies.\n", + "The bodies move separately in heave but move together in surge and pitch.\n", + "Therefore, we should define the individual heave inertia values _for each body_, but define the _total_ surge and pitch inertia values.\n", + "The inertia then needs to be reformatted as an `xarray.DataArray` to work with Capytaine.\n", + "The hydrostatic stiffness can be calculated for the total immersed body." ] }, { @@ -343,8 +370,23 @@ "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.\n", + "The solver will be run once for each wave phase realization.\n", + "Each of these phase realizations leads to a slightly different result for optimal average power.\n", + "Thus, for irregular wave conditions, it is recommended to include multiple phase realizations.\n", + "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", + "For this tutorial, the number of realizations has been set to 2 to reduce the total runtime. \n", + "\n", + "$$ t_{total} = \\frac{n_{realizations}}{f_1} $$\n", + "\n", + "Because we are now using irregular waves, we need significantly more frequencies to capture the entire wave spectrum and WEC response. " ] }, { @@ -353,11 +395,60 @@ "metadata": {}, "outputs": [], "source": [ - "# compute hydrodynamic coefficients\n", + "waves = {}\n", + "\n", "f1 = 0.02\n", "nfreq = 50\n", "freq = wot.frequency(f1, nfreq, False)\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 coefficients of the device, as done in previous tutorials.\n", + "Capytaine can handle generalized modes and will calculate the coefficients for our 4 degrees of freedom.\n", + "The BEM coefficients have been pre-calculated and are saved in a file.\n", + "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": [ "# read BEM data file if it exists\n", "filename = 'data/bem.nc'\n", "try:\n", @@ -371,7 +462,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now visualize the BEM results. An irregular frequency at about 5.25 rad/s is visible in many of the plots, but this should not substantially impact our model as the wave energy at this frequency is fairly low. " + "We now visualize the BEM results.\n", + "An irregular frequency at about 5.25 rad/s is visible in many of the plots, but this should not substantially impact our model as the wave energy at this frequency is fairly low. " ] }, { @@ -380,65 +472,7 @@ "metadata": {}, "outputs": [], "source": [ - "# plot coefficients\n", - "radiating_dofs = bem_data.radiating_dof.values\n", - "influenced_dofs = bem_data.influenced_dof.values\n", - "\n", - "# plots\n", - "fig_am, ax_am = plt.subplots(len(radiating_dofs), len(influenced_dofs),\n", - " tight_layout=True, sharex=True)\n", - "fig_rd, ax_rd = plt.subplots(len(radiating_dofs), len(influenced_dofs),\n", - " tight_layout=True, sharex=True)\n", - "fig_ex, ax_ex = plt.subplots(len(influenced_dofs), 1,\n", - " tight_layout=True, sharex=True, figsize=(4, 6))\n", - "\n", - "# plot titles\n", - "fig_am.suptitle('Added Mass Coefficients', fontweight='bold')\n", - "fig_rd.suptitle('Radiation Damping Coefficients', fontweight='bold')\n", - "fig_ex.suptitle('Wave Excitation Coefficients', fontweight='bold')\n", - "\n", - "# subplotting across 4DOF\n", - "sp_idx = 0\n", - "for i, rdof in enumerate(radiating_dofs):\n", - " for j, idof in enumerate(influenced_dofs):\n", - " sp_idx += 1\n", - " if i == 0:\n", - " np.abs(bem_data.diffraction_force.sel(influenced_dof=idof)).plot(\n", - " ax=ax_ex[j], linestyle='dashed', label='Diffraction force')\n", - " np.abs(bem_data.Froude_Krylov_force.sel(influenced_dof=idof)).plot(\n", - " ax=ax_ex[j], linestyle='dashdot', label='Froude-Krylov force')\n", - " ex_handles, ex_labels = ax_ex[j].get_legend_handles_labels()\n", - " ax_ex[j].set_title(f'{idof}')\n", - " ax_ex[j].set_xlabel('')\n", - " ax_ex[j].set_ylabel('')\n", - " if j <= i:\n", - " bem_data.added_mass.sel(\n", - " radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_am[i, j])\n", - " bem_data.radiation_damping.sel(\n", - " radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_rd[i, j])\n", - " if i == len(radiating_dofs)-1:\n", - " ax_am[i, j].set_xlabel(f'$\\omega$', fontsize=10)\n", - " ax_rd[i, j].set_xlabel(f'$\\omega$', fontsize=10)\n", - " ax_ex[j].set_xlabel(f'$\\omega$', fontsize=10)\n", - " else:\n", - " ax_am[i, j].set_xlabel('')\n", - " ax_rd[i, j].set_xlabel('')\n", - " if j == 0:\n", - " ax_am[i, j].set_ylabel(f'{rdof}', fontsize=10)\n", - " ax_rd[i, j].set_ylabel(f'{rdof}', fontsize=10)\n", - " else:\n", - " ax_am[i, j].set_ylabel('')\n", - " ax_rd[i, j].set_ylabel('')\n", - " if j == i:\n", - " ax_am[i, j].set_title(f'{idof}', fontsize=10)\n", - " ax_rd[i, j].set_title(f'{idof}', fontsize=10)\n", - " else:\n", - " ax_am[i, j].set_title('')\n", - " ax_rd[i, j].set_title('')\n", - " else:\n", - " fig_am.delaxes(ax_am[i, j])\n", - " fig_rd.delaxes(ax_rd[i, j])\n", - "fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False)" + "wot.utilities.plot_hydrodynamic_coefficients(bem_data)" ] }, { @@ -450,13 +484,13 @@ "The PTO model is similar to the one developed in Tutorial 2 but using the values corresponding to the LUPA PTO. \n", "The main difference is that in the LUPA the gear ratio can be modified by changing the interchangeable sprocket for one with a different radius. \n", "The motivation here is that operation in different wave conditions or different control schemes might have different torque and speed requirements. \n", - "Oregon State University has three sprockets of diameters 81.5mm (8MX-32S-36), 127.3mm (8MX-50S-36), 203.7mm (8MX-80S-36), but the manufacturer provides a fine selection of radius size in this range. \n", + "Oregon State University has three sprockets of diameters 81.5mm (8MX-32S-36), 127.3mm (8MX-50S-36), 203.7mm (8MX-80S-36), but the manufacturer provides a larger selection of radius size in this range. \n", "\n", "The PTO system consists of a [generator](https://akribis-systems.s3-us-west-2.amazonaws.com/pdfs/catalogs/adr-b.pdf), the [interchangeable sprocket](https://assets.gates.com/content/dam/gates/home/knowledge-center/resource-library/catalogs/old-pc_carbon_manual17595_2011.pdf), and two [idler pulleys](https://dpk3n3gg92jwt.cloudfront.net/domains/gates.pt/pdf/77234200.pdf) driven by a belt.\n", "\n", "

\"Schematic \"Digital

\n", "\n", - "We start by defining all the manufacturer-specified components." + "We start by defining all the manufacturer-specified components:" ] }, { @@ -617,7 +651,7 @@ "}\n", "\n", "# generator\n", - "# Note: Model ADR220-B175 data not listed online, but we have\n", + "# Note: Model ADR220-B175 data is not listed online, but we have a\n", "# spec sheet available on request \n", "generator = {\n", " 'torque_constant': 8.51, # N*m/A\n", @@ -702,9 +736,14 @@ "### Constraints\n", "In this tutorial we include realistic constraints on both the device motions and the generator operation:\n", "\n", - "* The **maximum stroke** (difference in heave between the float and spar) is 0.5m. Although there are end-stops (hard stop), the controller should ensure a soft stop. So, instead of modeling the hard stop we will add a constraint. \n", - "* The generator has both a **maximum torque** and **maximum speed**. It is important that we do not exceed that during operation.\n", - "* Generators also have a **continuous torque** requirement. The RMS torque during operation should not exceed this value or we risk damage to and failure of the machine. The torque can go higher that this (up to the maximum) for brief periods as long as the RMS torque remains below this threshold. \n", + "* The **maximum stroke** (difference in heave between the float and spar) is 0.5m.\n", + "Although there are end-stops (hard stop), the controller should ensure a soft stop.\n", + "So, instead of modeling the hard stop we will add a constraint. \n", + "* The generator has both a **maximum torque** and **maximum speed**.\n", + "It is important that we do not exceed that during operation.\n", + "* Generators also have a **continuous torque** requirement.\n", + "The RMS torque during operation should not exceed this value or we risk damage to and failure of the machine.\n", + "The torque can go higher that this (up to the maximum) for brief periods as long as the RMS torque remains below this threshold. \n", "\n", "We enforce all these through constraints on our design optimization problem. " ] @@ -759,7 +798,10 @@ "metadata": {}, "source": [ "### Mooring system\n", - "To fully capture the dynamics acting on LUPA in the Large Wave Flume, we must model the mooring system being used to account for the restoring forces acting on the WEC. The LUPA setup uses a 4-line taut mooring system with springs connecting the spar to the wall of the wave flume. The following mooring system properties have been provided by Oregon State University. The initial fairlead coordinates and anchor coordinates are relative to the center of gravity of the combined device:" + "To fully capture the dynamics acting on LUPA in the Large Wave Flume, we must model the mooring system being used to account for the restoring forces acting on the WEC.\n", + "The LUPA setup uses a 4-line taut mooring system with springs connecting the spar to the wall of the wave flume.\n", + "The following mooring system properties have been provided by Oregon State University.\n", + "The initial fairlead coordinates and anchor coordinates are relative to the center of gravity of the combined device:" ] }, { @@ -784,18 +826,31 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "There are several analytical and numerical methods commonly used to model mooring system kinematics for offshore systems, ranging from static analysis to determine equilibrium forces, all the way to FEA used to calculate the fully dynamic response of the mooring system components. For this design problem, we're mostly concerned with modeling the correct response of LUPA due to operational waves, so the high-fidelity methods are unnecessary at this design stage. While a purely linearized approach is common here, the symmetry of the taut lines in this current system allows us to instead use an analytical solution derived by Al-Solihat and Nahon (https://doi.org/10.1080/17445302.2015.1089052), which allows us to capture nonlinear mooring effects and off-diagonal terms in the mooring stiffness matrix without any significant increase in computation time.\n", + "There are several analytical and numerical methods commonly used to model mooring system kinematics for offshore systems, ranging from static analysis to determine equilibrium forces, all the way to FEA used to calculate the fully dynamic response of the mooring system components.\n", + "For this design problem, we are mostly concerned with modeling the correct response of LUPA due to operational waves, so the high-fidelity methods are unnecessary at this design stage.\n", + "While a purely linearized approach is common here, the symmetry of the taut lines in this current system allows us to instead use an analytical solution derived by Al-Solihat and Nahon ([link](https://doi.org/10.1080/17445302.2015.1089052)), which allows us to capture nonlinear mooring effects and off-diagonal terms in the mooring stiffness matrix without any significant increase in computation time.\n", "\n", "This solution takes an exact analysis of the derivatives of the classic elastic catenary equations and simplifies them by assuming the taut lines have no sag and negligible mass, allowing for the differential changes of the horizontal and vertical restoring force to be calculated as\n", + "\n", "$$ \\frac{\\partial F_{hor}}{\\partial l} = K_{axial}\\textup{cos}^2 \\theta + \\frac{T}{L} \\textup{sin}^2 \\theta $$\n", "$$ \\frac{\\partial F_{vert}}{\\partial h} = K_{axial}\\textup{sin}^2 \\theta + \\frac{T}{L} \\textup{cos}^2 \\theta $$\n", "$$ \\frac{\\partial F_{hor}}{\\partial h} = \\frac{\\partial F_{vert}}{\\partial l} = \\textup{cos} \\theta \\textup{sin} \\theta \\left [ K_{axial}- \\frac{T}{L}\\right ]$$\n", "\n", - "where $l$ and $h$ are the horizontal and vertical distance between the anchor points and fairlead points, respectively; $T$ is the pretension; $\\theta$ is the angle between the seabed and the mooring line such that $ \\theta = \\textup{tan}^{-1}(\\frac{h}{l}) $; and $L$ is the stretched length of the mooring line such that $L = \\sqrt{l^2 + h^2} $\n", + "where $l$ and $h$ are the horizontal and vertical distance between the anchor points and fairlead points, respectively; $T$ is the pretension; $\\theta$ is the angle between the seabed and the mooring line such that\n", + "\n", + "$$ \\theta = \\textup{tan}^{-1}(\\frac{h}{l}) $$\n", + "\n", + "and $L$ is the stretched length of the mooring line such that\n", + "\n", + "$$ L = \\sqrt{l^2 + h^2} $$\n", "\n", "When these equations are applied to the linear stiffness equation in each radiating ($i$) and influencing ($j$) degrees of freedom:\n", - "$$\\boldsymbol{K_{mooring}} = - \\frac{\\partial \\boldsymbol{F_{mooring}}}{\\partial \\boldsymbol{X}} = \\sum_{m=1}^{n_{lines}} [K_{ij}]^{(m)} = - \\sum_{m=1}^{n_{lines}}[\\frac{\\partial (F_{mooring})_i}{\\partial X_j}]^{(m)} $$\n", - "where $X$ are the generalized displacements of the WEC in each degree of freedom, they yield Equation (27) from the reference above which provides an analytical solution to the mooring stiffness matrix, which translates to the `k_mooring` function below. See the reference above for the full theoretical explanation and derivation of these equations." + "\n", + "$$\\boldsymbol{K_{mooring}} = - \\frac{\\partial \\boldsymbol{F_{mooring}}}{\\partial \\boldsymbol{X}} \\\\\n", + "= \\sum_{m=1}^{n_{lines}} [K_{ij}]^{(m)} = - \\sum_{m=1}^{n_{lines}}[\\frac{\\partial (F_{mooring})_i}{\\partial X_j}]^{(m)} $$\n", + "\n", + "where $X$ are the generalized displacements of the WEC in each degree of freedom, they yield Equation (27) from the reference above which provides an analytical solution to the mooring stiffness matrix, which translates to the `k_mooring` function below.\n", + "See the reference above for the full theoretical explanation and derivation of these equations." ] }, { @@ -855,7 +910,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "`k_mooring` defines the 7x7 mooring matrix for a two-body point absorber WEC. Since we are only concerned with the planar components here, we can extract the rows/columns for the heaves, surge, and pitch to obtain the 4x4 matrix we need here." + "`k_mooring` defines the 7x7 mooring matrix for a two-body point absorber WEC.\n", + "Since we are only concerned with the planar components here, we can extract the rows/columns for the heaves, surge, and pitch to obtain the 4x4 matrix we need here." ] }, { @@ -893,7 +949,7 @@ "tmp = moor.isel(omega=0).copy(deep=True)\n", "tmp['omega'] = tmp['omega'] * 0\n", "moor = xr.concat([tmp, moor], dim='omega') \n", - "moor = moor.transpose(\"omega\", \"radiating_dof\", \"influenced_dof\")\n", + "moor = moor.transpose('omega', 'radiating_dof', 'influenced_dof')\n", "moor = -1*moor # RHS of equation: -ma = ÎŁf \n", "mooring_force = wot.force_from_rao_transfer_function(moor, True)" ] @@ -902,16 +958,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Pre-tension of spar\n", + "#### Pretension of spar\n", "\n", - "The spar is positively buoyant at its equilibrium position and relies on pre-tension from the mooring. \n", - "The mooring matrix above only captures the restoring forces from the mooring lines, not the pretension itself. \n", - "Here we have two choices:\n", - "\n", - "1. We can ignore this fact. WecOptTool asssumes equilibrium so the pre-tension force is implied. \n", - "2. Explicitly model the buoyancy, gravity, and pre-tension forces as done in Tutorial 2. \n", - "\n", - "In Tutorial 2 we had constraints on the mooring/tether line (minimum tension, never slack) and therefore it was important to explicitly model these forces. In this case we do not have such constraints and will simply allow to pre-tension force to be implied. This does not affect the solution." + "The spar is positively buoyant at its equilibrium position and relies on pretension from the mooring. \n", + "The mooring matrix above only captures the restoring forces from the mooring lines, not the pretension itself.\n", + "However, unlike Tutorial 2, we can ignore modeling this; WecOptTool assumes the device is starting in equilibrium, thus the pretension force is implied.\n", + "We modeled this in Tutorial 2 because that optimization problem (comparing the hull mass vs. pretension) required relevant constraints on the mooring/tether line, so we had to explicitly model the buoyancy, gravity, and pretension force there to make those constraints possible.\n", + "These forces have no impact on the solution to the optimization problem of interest in Part 2 of this tutorial (optimal sprocket sizing), so we can ignore them here." ] }, { @@ -947,101 +1000,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": {}, - "source": [ - "We now plot the wave spectra in order to verify that our frequency array is sufficient to capture our target waves. The peak frequency for the PacWave South maximum occurrence is cut off somewhat, so we would want to sample more wave frequencies (i.e. increase `nfreq`) in a formal analysis to properly capture this. However, this is sufficient for purposes of this tutorial." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots()\n", - "plt1 = np.abs(waves['south_max_90'].sel(realization=0)).plot(\n", - " ax=ax, color='C0', linestyle='solid', label='PW South, 90th percentile')\n", - "plt2 = np.abs(waves['south_max_annual'].sel(realization=0)).plot(\n", - " ax=ax, color='C0', linestyle='dotted', label='PW South, Max Annual')\n", - "plt3 = np.abs(waves['south_max_occurrence'].sel(realization=0)).plot(\n", - " ax=ax, color='C0', linestyle='dashed', label='PW South, Max Occurrence')\n", - "plt4 = np.abs(waves['south_min_10'].sel(realization=0)).plot(\n", - " ax=ax, color='C0', linestyle='dashdot', label='PW South, 10th percentile')\n", - "plt5 = np.abs(waves['north_max_90'].sel(realization=0)).plot(\n", - " ax=ax, color='C1', linestyle='solid', label='PW North, 90th percentile')\n", - "plt6 = np.abs(waves['north_max_annual'].sel(realization=0)).plot(\n", - " ax=ax, color='C1', linestyle='dotted', label='PW North, Max Annual')\n", - "plt7 = np.abs(waves['north_max_occurrence'].sel(realization=0)).plot(\n", - " ax=ax, color='C1', linestyle='dashed', label='PW North, Max Occurrence')\n", - "plt8 = np.abs(waves['north_min_10'].sel(realization=0)).plot(\n", - " ax=ax, color='C1', linestyle='dashdot', label='PW North, 10th percentile')\n", - "\n", - "ax.set_title('PacWave wave spectra, LWF scale', fontweight='bold')\n", - "plts = plt1 + plt2 + plt3 + plt4 + plt5 + plt6 + plt7 + plt8\n", - "ax.legend(plts, [pl.get_label() for pl in plts], ncols=1, frameon=False)" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -1049,7 +1007,9 @@ "### Solve\n", "We now solve the inner optimization for optimal control strategy for a fixed design.\n", "\n", - "In a formal design case, it is good practice to first test the setup with a regular wave case to confirm the WEC is responding as expected. We have performed this step, but skip it here for the sake of brevity. To run the regular wave, change the wave selection in the code below. " + "In a formal design case, it is good practice to first test the setup with a regular wave case to confirm the WEC is responding as expected.\n", + "We have performed this step, but skip it here for the sake of brevity.\n", + "To run the regular wave, change the wave selection in the code below. " ] }, { @@ -1068,8 +1028,7 @@ "scale_obj = 1e-2\n", "\n", "results = wec.solve(\n", - " waves['south_max_occurrence'],\n", - " # waves['regular'], \n", + " waves['south_max_occurrence'], #waves['regular'],\n", " obj_fun, \n", " nstate_opt, \n", " scale_x_wec=scale_x_wec,\n", @@ -1080,16 +1039,6 @@ "print(f'Optimal average electrical power: {np.mean(power_results)} W')" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x_wec, x_opt = wec.decompose_state(results[0].x)\n", - "print(len(x_opt))" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1098,8 +1047,8 @@ "source": [ "# Post-process\n", "nsubsteps = 5\n", - "wec_fdom, wec_tdom = wec.post_process(results[0], waves['south_max_occurrence'].sel(realization=0), nsubsteps)\n", - "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves['south_max_occurrence'].sel(realization=0), nsubsteps)" + "wec_fdom, wec_tdom = wec.post_process(wec, results, waves['south_max_occurrence'], nsubsteps)\n", + "pto_fdom, pto_tdom = pto.post_process(wec, results, waves['south_max_occurrence'], nsubsteps)" ] }, { @@ -1107,11 +1056,17 @@ "metadata": {}, "source": [ "### Results\n", - "Looking at the response of the optimal solution, we can see the behavior of the WEC and its relationship to our problem constraints:\n", - "\n", - "* __PTO position__ - The PTO position corresponds to the difference in the vertical position between the two bodies. The PTO position is clearly nowhere near the 0.5 maximum stroke, so the constraint is not plotted here.\n", - "* __Excitation force and velocity__ - Much like Tutorial 2, we would expect the excitation force and velocity of each body to be in phase when maximizing for mechanical power and in the absence of constraints. Given the constraints on the PTO, this is not the case here, with the spar velocity being out of phase with the wave excitation acting on it. This indicates the differences when maximizing for mechanical vs. electrical power for LUPA in these wave conditions. It is also clear that the velocities (or their rotational equivalents) are nowhere near the generator maximum speed.\n", - "* __Generator torque__ - This plot shows that both of the torque constraints are active for this wave case, and dominate the solution given the lack of relevance of the other constraints above. The PTO RMS torque is well within the continuous torque limit, and the generator exceeds the RMS limit for only brief period. The peak torque limit is approached at a few moments as well.\n", + "The `post_process` functions return lists of post-processed results. Here, we extract the first element of each list to analyze the results for the first wave phase realization. Looking at the response of the optimal solution, we can see the behavior of the WEC and its relationship to our problem constraints:\n", + "\n", + "* __PTO position__ - The PTO position corresponds to the difference in the vertical position between the two bodies.\n", + "The PTO position is clearly nowhere near the 0.5 m maximum stroke, so the constraint is not plotted here.\n", + "* __Excitation force and velocity__ - Much like Tutorial 2, we would expect the excitation force and velocity of each body to be in phase when maximizing for mechanical power and in the absence of constraints.\n", + "Given the constraints on the PTO, this is not the case here, with the spar velocity being out of phase with the wave excitation acting on it.\n", + "This indicates the differences when maximizing for mechanical vs. electrical power for LUPA in these wave conditions.\n", + "It is also clear that the velocities (or their rotational equivalents) are nowhere near the generator maximum speed.\n", + "* __Generator torque__ - This plot shows that both of the torque constraints are active for this wave case, and dominate the solution given the lack of relevance of the other constraints above.\n", + "The PTO RMS torque is well within the continuous torque limit, and the generator exceeds the RMS limit for only brief period.\n", + "The peak torque limit is approached at a few moments as well.\n", "* __Power__ - As expected, the greatest power absorption aligns with the times of greatest PTO velocity and peak torque magnitude. " ] }, @@ -1124,11 +1079,11 @@ "fig_res1, ax_res1 = plt.subplots(nrows=4, sharex=True, figsize=(8, 10))\n", "\n", "## PTO position\n", - "wec_tdom.pos.sel(influenced_dof='spar__Heave').plot(\n", + "wec_tdom[0].pos.sel(influenced_dof='spar__Heave').plot(\n", " ax=ax_res1[0], linestyle='dashed', color='C7', label='Spar')\n", - "wec_tdom.pos.sel(influenced_dof='float__Heave').plot(\n", + "wec_tdom[0].pos.sel(influenced_dof='float__Heave').plot(\n", " ax=ax_res1[0], linestyle='dotted', color='C6', label='Float')\n", - "pto_tdom.pos.plot(ax=ax_res1[0], label='PTO', color='black')\n", + "pto_tdom[0].pos.plot(ax=ax_res1[0], label='PTO', color='black')\n", "\n", "ax_res1[0].set_ylabel('Position [m]')\n", "ax_res1[0].set_title('PTO position vs. WEC body position', fontweight='bold')\n", @@ -1141,23 +1096,23 @@ "twinax = ax_res1[1].twinx()\n", "\n", "# Spar excitation\n", - "force_excitation_spar = wec_tdom.force.sel(\n", + "force_excitation_spar = wec_tdom[0].force.sel(\n", " influenced_dof='spar__Heave', type=['Froude_Krylov', 'diffraction'])\n", "force_excitation_spar = force_excitation_spar.sum('type')\n", "plt1 = force_excitation_spar.plot(\n", " ax=ax_res1[1], linestyle='dashed', color='C0', label='Spar exc.')\n", "\n", "# Float excitation\n", - "force_excitation_float = wec_tdom.force.sel(\n", + "force_excitation_float = wec_tdom[0].force.sel(\n", " influenced_dof='float__Heave', type=['Froude_Krylov', 'diffraction'])\n", "force_excitation_float = force_excitation_float.sum('type')\n", "plt2 = force_excitation_float.plot(\n", " ax=ax_res1[1], linestyle='dotted', color='C0', label='Float exc.')\n", "\n", "# Spar/float velocity\n", - "plt3 = wec_tdom.vel.sel(influenced_dof='spar__Heave').plot(\n", + "plt3 = wec_tdom[0].vel.sel(influenced_dof='spar__Heave').plot(\n", " ax=twinax, color='C1', linestyle='dashed', label='Spar velocity')\n", - "plt4 = wec_tdom.vel.sel(influenced_dof='float__Heave').plot(\n", + "plt4 = wec_tdom[0].vel.sel(influenced_dof='float__Heave').plot(\n", " ax=twinax, color='C1', linestyle='dotted', label='Float velocity')\n", "\n", "twinax.set_ylabel('Velocity [m/s]', color='C1')\n", @@ -1175,26 +1130,26 @@ "twinax.set_ylim(bottom=twinax.get_ylim()[0]*1.25)\n", "\n", "## Torque\n", - "pto_tdom.force.plot(\n", + "pto_tdom[0].force.plot(\n", " ax=ax_res1[2], linestyle='solid', color='black', label='PTO trq.')\n", - "pto_rms_tq = np.sqrt(np.mean(pto_tdom.force.values**2)\n", - " / pto_tdom.time.values[-1])\n", + "pto_rms_tq = np.sqrt(np.mean(pto_tdom[0].force.values**2)\n", + " / pto_tdom[0].time.values[-1])\n", "ax_res1[2].plot(\n", - " pto_tdom.time, pto_rms_tq * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, pto_rms_tq * np.ones(pto_tdom[0].time.shape),\n", " color='C9', linestyle='solid', label='PTO RMS trq.')\n", "max_tq = generator['max_torque']\n", "rms_tq = generator['continuous_torque']\n", "ax_res1[2].plot(\n", - " pto_tdom.time, 1*max_tq * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, 1*max_tq * np.ones(pto_tdom[0].time.shape),\n", " color='C5', linestyle='dotted', label=f'Peak trq. limit')\n", "ax_res1[2].plot(\n", - " pto_tdom.time, -1*max_tq * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, -1*max_tq * np.ones(pto_tdom[0].time.shape),\n", " color='C5', linestyle='dotted')\n", "ax_res1[2].plot(\n", - " pto_tdom.time, 1*rms_tq * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, 1*rms_tq * np.ones(pto_tdom[0].time.shape),\n", " color='C9', linestyle='dotted', label=f'RMS trq. limit')\n", "ax_res1[2].plot(\n", - " pto_tdom.time, -1*rms_tq * np.ones(pto_tdom.time.shape),\n", + " pto_tdom[0].time, -1*rms_tq * np.ones(pto_tdom[0].time.shape),\n", " color='C9', linestyle='dotted')\n", "ax_res1[2].set_ylabel('Torque [Nm] ')\n", "ax_res1[2].legend(ncols=4, loc='lower center', frameon=False)\n", @@ -1204,9 +1159,9 @@ "ax_res1[2].set_ylim(bottom=ax_res1[2].get_ylim()[0]*1.5)\n", "\n", "## Power\n", - "(pto_tdom['power'].loc['mech',:,:]).plot(\n", + "(pto_tdom[0]['power'].loc['mech',:,:]).plot(\n", " ax=ax_res1[3], color='C8', label='Mech. power')\n", - "(pto_tdom['power'].loc['elec',:,:]).plot(\n", + "(pto_tdom[0]['power'].loc['elec',:,:]).plot(\n", " ax=ax_res1[3], color='C5', linestyle='dashed', label=\"Elec. power\")\n", "ax_res1[3].legend(ncols=2, loc='lower center', frameon=False)\n", "ax_res1[3].set_title('Power generation', fontweight='bold')\n", @@ -1219,7 +1174,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can also look at all four degrees of freedom and make sure the relationships between the forces in each direction are reasonable. We can see that the PTO and radiation forces are generally in phase in the heave direction, but inversely related for the spar. We also see that the mooring response is inversely propotional to the radiation forces of the spar (with a more delayed restoring response in pitch than in surge/heave), and that the mooring correctly has no effect on the float, which it doesn't connect to." + "We can also look at all four degrees of freedom and make sure the relationships between the forces in each direction are reasonable.\n", + "We can see that the PTO and radiation forces are generally in phase in the heave direction, but inversely related for the spar.\n", + "We also see that the mooring response is inversely propotional to the radiation forces of the spar (with a more delayed restoring response in pitch than in surge/heave), and that the mooring correctly has no effect on the float, which it does not connect to." ] }, { @@ -1228,21 +1185,22 @@ "metadata": {}, "outputs": [], "source": [ + "influenced_dofs = bem_data.influenced_dof.values\n", "## Forces\n", "fig_res2, ax_res2 = plt.subplots(len(influenced_dofs), 1, sharex=True, figsize=(8, 10))\n", "for j, idof in enumerate(influenced_dofs):\n", - " force_excitation = wec_tdom.force.sel(\n", + " force_excitation = wec_tdom[0].force.sel(\n", " type=['Froude_Krylov', 'diffraction'])\n", " force_excitation = force_excitation.sum('type')\n", " force_excitation.sel(influenced_dof=idof).plot(\n", " ax=ax_res2[j], label='Excitation')\n", - " wec_tdom.force.sel(influenced_dof=idof, type='radiation').plot(\n", + " wec_tdom[0].force.sel(influenced_dof=idof, type='radiation').plot(\n", " ax=ax_res2[j], label='Radiation')\n", - " wec_tdom.force.sel(influenced_dof=idof, type='hydrostatics').plot(\n", + " wec_tdom[0].force.sel(influenced_dof=idof, type='hydrostatics').plot(\n", " ax=ax_res2[j], label='Hydrostatics')\n", - " wec_tdom.force.sel(influenced_dof=idof, type='Mooring').plot(\n", + " wec_tdom[0].force.sel(influenced_dof=idof, type='Mooring').plot(\n", " ax=ax_res2[j], label='Mooring')\n", - " wec_tdom.force.sel(influenced_dof=idof, type='PTO').plot(\n", + " wec_tdom[0].force.sel(influenced_dof=idof, type='PTO').plot(\n", " ax=ax_res2[j], label='PTO')\n", " ax_res2[j].set_title(f'{idof}')\n", " if j != len(influenced_dofs)-1:\n", @@ -1257,15 +1215,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Control co-design of the PTO sprocket sizing for maximum electrical power" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ + "## 2. Control co-design of the PTO sprocket sizing for maximum electrical power\n", + "\n", "### Setup\n", - "With our model working as expected, we can now iterate on this model with varying sprocket sizing to identify the optimal size. As with previous tutorials, we wrap the code from Part 1 into a function and iterate on our chosen design parameter (i.e. each key of the `sprockets` dictionary earlier corresponds to the `x` argument of `design_obj_fun`)" + "With our model working as expected, we can now iterate on this model with varying sprocket sizing to identify the optimal size.\n", + "As with previous tutorials, we wrap the code from Part 1 into a function and iterate on our chosen design parameter (i.e. each key of the `sprockets` dictionary earlier corresponds to the `x` argument of `design_obj_fun`):" ] }, { @@ -1350,8 +1304,7 @@ " f' Sprocket design: {spr_des}')\n", "\n", " results = wec.solve(\n", - " waves['south_max_occurrence'],\n", - " # waves['regular'], \n", + " waves['south_max_occurrence'], #waves['regular'],\n", " obj_fun, \n", " nstate_opt, \n", " scale_x_wec=scale_x_wec,\n", @@ -1379,7 +1332,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Given there are only 22 sprockets, we will continue to use a brute force algorithm. This algorithm typically takes several hours to run; for sake of time here, we have included the calculated results in the `tutorial_3_results.nc` file. If you would like to run the brute force algorithm yourself (e.g. if you modify this notebook and would like to see how the results change), move or delete this file from the `data` directory." + "Given there are only 22 sprockets, we will continue to use a brute force algorithm.\n", + "This algorithm typically takes several hours to run; for sake of time here, we have included the calculated results in the `tutorial_3_results.nc` file.\n", + "If you would like to run the brute force algorithm yourself (e.g. if you modify this notebook and would like to see how the results change), move or delete this file from the `data` directory." ] }, { @@ -1416,9 +1371,13 @@ "metadata": {}, "source": [ "### Results\n", - "The power generation for all the sprockets are plotted below, with the three sprockets possessed by Oregon State University called out. Comparing the results across the range of sprocket sizes, it is clear that power generation is maximized towards smaller sprocket diameters. Here, 8MX-33S-36 sprocket (wiht the second smallest diameter) is the best sprocket selection for LUPA for the selected wave case.\n", + "The power generation for all the sprockets are plotted below, with the three sprockets possessed by Oregon State University called out.\n", + "Comparing the results across the range of sprocket sizes, it is clear that power generation is maximized towards smaller sprocket diameters.\n", + "Here, 8MX-33S-36 sprocket (the second smallest diameter) is the best sprocket selection for LUPA for the selected wave case.\n", "\n", - "The motivation to include an interchangeable sprocket in LUPA is to make it versatile across different locations and wave conditions. The optimal sprocket could be different when tested with a more or less severe wave climate, or a different mooring configuration. Try moving the `data/tutorial_3_results.nc` file as above, and test different `wave_cases` keys or modify the mooring system parameters (`init_fair_coords`, `anch_coords`, `line_ax_stiff`, and `pretension`) and see how the WEC response changes and which sprocket generates the most power." + "The motivation to include an interchangeable sprocket in LUPA is to make it versatile across different locations and wave conditions.\n", + "The optimal sprocket could be different when tested with a more or less severe wave climate, or a different mooring configuration.\n", + "Try moving the `data/tutorial_3_results.nc` file as above, and test different `wave_cases` keys or modify the mooring system parameters (`init_fair_coords`, `anch_coords`, `line_ax_stiff`, and `pretension`) and see how the WEC response changes and which sprocket generates the most power." ] }, { diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb index a817b7a7b..d7aa24dcd 100644 --- a/examples/tutorial_4_Pioneer.ipynb +++ b/examples/tutorial_4_Pioneer.ipynb @@ -5,22 +5,27 @@ "metadata": {}, "source": [ "# Tutorial 4 - Pioneer\n", - "This tutorial models a potential design for a WEC used to provide power to instruments in the [Pioneer Central Surface Mooring System](https://oceanobservatories.org/site/cp01cnsm/) within the [Coastal Pioneer Array](https://oceanobservatories.org/array/coastal-pioneer-array/) off the coast of New England. This system has instrumentation to measure various meteorological, surface, near-surface, and seabed phenomena on the Continental Shelf-Slope and transmit back to shore.\n", + "This tutorial models a potential design for a WEC used to provide power to instruments in the [Pioneer Central Surface Mooring System](https://oceanobservatories.org/site/cp01cnsm/) within the National Science Foundation Ocean Observatories Initiative [Coastal Pioneer Array](https://oceanobservatories.org/array/coastal-pioneer-array/).\n", + "This system has instrumentation to measure various meteorological, surface, near-surface, and seabed phenomena on the Continental Shelf-Slope and transmit back to shore.\n", "\n", - "Unlike previous tutorials, this tutorial does not include an outer optimization loop. Instead, the focus here is on effectively modeling an optimization problem for the control of a unique WEC archetype with a pitch resonator PTO with several co-dependent components. This is completely setup in Part 1, and includes:\n", + "Unlike previous tutorials, this tutorial does not include an outer optimization loop.\n", + "Instead, the focus here is on effectively modeling an optimization problem for the control of a unique WEC archetype with a pitch resonator PTO with several co-dependent components.\n", + "This is completely setup in Part 1, and includes:\n", "\n", "* Expanding the control state vector\n", "* Custom PTO physics\n", "* Modeling non-hydrodynamic degrees of freedom\n", "\n", - "This tutorial is divided into three parts. The first sets up the problem including the points above. The second and third show results using a regular and irregular wave, respectively.\n", + "This tutorial is divided into three parts.\n", + "The first sets up the problem including the points above.\n", + "The second and third show results using a regular and irregular wave, respectively.\n", "\n", "1. [Model setup](#1.-Model-setup)\n", "2. [Regular wave results](#2.-Regular-wave-results)\n", "3. [Irregular wave](#3.-Irregular-wave)\n", "\n", "\n", - "

\"Diagram

" + "

\"Diagram

" ] }, { @@ -33,7 +38,6 @@ "import autograd.numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.linalg import block_diag\n", - "import xarray as xr\n", "\n", "import wecopttool as wot\n", "\n", @@ -48,12 +52,15 @@ "## 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 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", + "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.\n", + "The regular wave is roughly at 0.35 Hz, the known pitch resonance frequency of the buoy.\n", + "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.\n", + "For the purpose of the tutorial, the number of realizations has been set to 2 to reduce runtime.\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 +69,9 @@ "metadata": {}, "outputs": [], "source": [ - "f1 = 0.025 # Hz\n", - "nfreq = 75" + "fend = 1.875\n", + "nfreq = 150\n", + "f1 = fend / nfreq" ] }, { @@ -81,24 +89,19 @@ "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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We plot the wave elevation spectra to check if the chosen frequency array covers the main excitation frequency, converges to zero for larger frequencies and can capture multiple harmonics which are required for nonlinear dynamics, since the states are multiplied with each other.\n", - "If the plot below does not cover the desired frequencies please adjust the frequency array." + "We plot the wave elevation spectra to check if the chosen frequency array covers the main excitation frequency, converges to zero for larger frequencies and can capture multiple harmonics which are required for nonlinear dynamics, since the states are multiplied with each other." ] }, { @@ -120,8 +123,10 @@ "metadata": {}, "source": [ "### 1.2 Basic properties\n", + "\n", "#### Geometry\n", - "To get the required hydrodynamic coefficients of the Pioneer design, we only need to be concerned with creating a mesh of the buoy geometry. The buoy has the same general shape as the WaveBot used in Tutorial 1, so we will use the build-in `geom.WaveBot` object and change the dimensions to simplify the mesh creation process." + "To get the required hydrodynamic coefficients of the Pioneer design, we only need to be concerned with creating a mesh of the buoy geometry.\n", + "The buoy has the same general shape as the WaveBot used in Tutorial 1, so we will use the build-in `geom.WaveBot` object and change the dimensions to simplify the mesh creation process." ] }, { @@ -148,11 +153,18 @@ "metadata": {}, "source": [ "#### Design properties\n", - "The Pioneer team has developed initial specifications for a full-sized WEC design that fits within the buoy dimensions and operating at resonance. There are several components that we need to capture in the PTO model:\n", - "\n", - "* **Buoy**: As mentioned, this has the same basic geometric parameters as the WaveBot used in Tutorial 1. Since we're only modeling the pitch degree of freedom, we only need the pitch moment of inertia and not the mass of the buoy.\n", - "* **Flywheel**: The flywheel is contained and dimensioned to fit on top of the buoy. It rotates in pitch independently of the buoy, but is linked to the buoy through two parallel connections: the torsional spring and the generator. Since the flywheel is protected from wave loading by the buoy, it is **not** subject to hydrostatics or hydrodynamics.\n", - "* **Spring**: The magnetic torsional spring provides positive stiffness to the flywheel. It has a linear displacement up to 45 degrees, at which point it \"wraps around\" and the stiffness begins decreasing. For simplicity, we will model the spring as strictly linear. The spring is connected to a gearbox which makes the position of the spring (and therefore, the generator shaft) different from the relative position of the flywheel and buoy. \n", + "The Pioneer team has developed initial specifications for a full-sized WEC design that fits within the buoy dimensions and operating at resonance.\n", + "There are several components that we need to capture in the PTO model:\n", + "\n", + "* **Buoy**: As mentioned, this has the same basic geometric parameters as the WaveBot used in Tutorial 1.\n", + "Since we're only modeling the pitch degree of freedom, we only need the pitch moment of inertia and not the mass of the buoy.\n", + "* **Flywheel**: The flywheel is contained and dimensioned to fit on top of the buoy.\n", + "It rotates in pitch independently of the buoy, but is linked to the buoy through two parallel connections: the torsional spring and the generator.\n", + "Since the flywheel is protected from wave loading by the buoy, it is **not** subject to hydrostatics or hydrodynamics.\n", + "* **Spring**: The magnetic torsional spring provides positive stiffness to the flywheel.\n", + "It has a linear displacement up to 45 degrees, at which point it \"wraps around\" and the stiffness begins decreasing.\n", + "For simplicity, we will model the spring as strictly linear.\n", + "The spring is connected to a gearbox which makes the position of the spring (and therefore, the generator shaft) different from the relative position of the flywheel and buoy. \n", "\n", "A rendering of the Pioneer WEC model is shown here with the buoy pitch in red and the flywheel pitch in green:\n", "

\"Diagram

" @@ -188,15 +200,22 @@ "metadata": {}, "source": [ "#### Pitch resonator\n", - "This system uses a unique PTO system, where the generator velocity is dependent on the relative pitching rotation of the buoy and the flywheel. The system can be modeled as a coupled spring-mass-damper system, with the generator represented as $B_{PTO}$. The spring-mass-damper system and equivalent circuit model are shown below. Note that all of these values are rotational and are drawn here linearly for simplicity. \n", + "This system uses a unique PTO system, where the generator velocity is dependent on the relative pitching rotation of the buoy and the flywheel.\n", + "The system can be modeled as a coupled spring-mass-damper system, with the generator represented as $B_{PTO}$.\n", + "The spring-mass-damper system and equivalent circuit model are shown below.\n", + "Note that all of these values are rotational and are drawn here linearly for simplicity. \n", "\n", "

\"Pioneer \"Pioneer

\n", "\n", "This corresponds to the following dynamics equations:\n", "\n", - "(1) $$ [J_{buoy} + J(\\omega))] \\ddot{\\theta}_{buoy} = T_{hydro} - K_{hydro}\\theta_{buoy} - B(\\omega)\\dot{\\theta}_{buoy} - K_{spring}(\\theta_{buoy} - \\theta_{flywheel})/N^2 - (B_{PTO} + B_{fric})(\\dot{\\theta}_{buoy} - \\dot{\\theta}_{flywheel}) $$\n", + "(1)\n", + "\n", + "$$ [J_{buoy} + J(\\omega))] \\ddot{\\theta}_{buoy} = T_{hydro} - K_{hydro}\\theta_{buoy} - B(\\omega)\\dot{\\theta}_{buoy} - K_{spring}(\\theta_{buoy} - \\theta_{flywheel})/N^2 - (B_{PTO} + B_{fric})(\\dot{\\theta}_{buoy} - \\dot{\\theta}_{flywheel}) $$\n", "\n", - "(2) $$ J_{flywheel} \\ddot{\\theta}_{flywheel} = K_{spring}(\\theta_{buoy} - \\theta_{flywheel})/N^2 + (B_{PTO} + B_{fric})(\\dot{\\theta}_{buoy} - \\dot{\\theta}_{flywheel}) $$\n", + "(2)\n", + "\n", + "$$ J_{flywheel} \\ddot{\\theta}_{flywheel} = K_{spring}(\\theta_{buoy} - \\theta_{flywheel})/N^2 + (B_{PTO} + B_{fric})(\\dot{\\theta}_{buoy} - \\dot{\\theta}_{flywheel}) $$\n", "\n", "We will explore how to capture these dynamics in WecOptTool in the following sections." ] @@ -206,7 +225,7 @@ "metadata": {}, "source": [ "#### Hydrodynamics and hydrostatics\n", - "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves. The inertia matrix needs to be defined manually here since it is not based on the displaced mass." + "As mentioned above, the `FloatingBody` object in Capytaine only needs to model the buoy, since no other components are being excited by the waves." ] }, { @@ -220,32 +239,24 @@ "pnr_fb.center_of_mass = np.array([0., 0., buoy_props['CG']])\n", "pnr_fb.rotation_center = pnr_fb.center_of_mass\n", "ndof = pnr_fb.nb_dofs\n", - "pnr_fb.show_matplotlib()\n", - "\n", - "pnr_fb.inertia_matrix = xr.DataArray([[buoy_props['MOI']]],\n", - " dims=['influenced_dof', 'radiating_dof'],\n", - " coords={'influenced_dof' : ['Pitch'],\n", - " 'radiating_dof' : ['Pitch']},\n", - " name=\"inertia_matrix\"\n", - " )" + "pnr_fb.show_matplotlib()" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "rho = 1025. # kg/m^3\n", "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n", - "bem_data = wot.run_bem(pnr_fb, freq, rho=rho)\n", + "bem_data = wot.run_bem(pnr_fb, freq)\n", "omega = bem_data.omega.values\n", "\n", "pnr_fb.keep_immersed_part()\n", + "k_buoy = pnr_fb.compute_hydrostatic_stiffness(rho=rho).values.squeeze()\n", "k_spring = spring_props['Max torque'] / spring_props['Max displacement']\n", - "print(f\"Hydrostatic stiffness from Capytaine: {bem_data['hydrostatic_stiffness'].values.squeeze()} N-m/rad\")\n", + "print(f'Hydrostatic stiffness from Capytaine: {k_buoy} N-m/rad')\n", "print('Hydrostatic stiffness from experiment: 37204 N-m/rad')" ] }, @@ -275,6 +286,7 @@ "ex_handles, ex_labels = ax_ex.get_legend_handles_labels()\n", "ax_ex.set_xlabel(f'$\\omega$', fontsize=10)\n", "ax_ex.set_title('Wave Excitation Coefficients', fontweight='bold')\n", + "ax_ex.grid(True)\n", "fig_ex.legend(ex_handles, ex_labels, loc='center right', frameon=False)\n", "\n", "# Added mass\n", @@ -282,12 +294,14 @@ " radiating_dof='Pitch', influenced_dof='Pitch').plot(ax=ax_am)\n", "ax_am.set_xlabel(f'$\\omega$', fontsize=10)\n", "ax_am.set_title('Added Mass Coefficients', fontweight='bold')\n", + "ax_am.grid(True)\n", "\n", "# Radiation damping\n", "bem_data.radiation_damping.sel(\n", " radiating_dof='Pitch', influenced_dof='Pitch').plot(ax=ax_rd)\n", "ax_rd.set_xlabel(f'$\\omega$', fontsize=10)\n", - "ax_rd.set_title('Radiation Damping Coefficients', fontweight='bold')" + "ax_rd.set_title('Radiation Damping Coefficients', fontweight='bold')\n", + "ax_rd.grid(True)" ] }, { @@ -295,10 +309,12 @@ "metadata": {}, "source": [ "### 1.3 PTO\n", - "The generator shaft velocity is the difference between the buoy velocity and the flywheel velocity. Because the velocity PTO depends not only on the hydrodynamic states (`x_wec`, buoy pitch) but also on an additional non-hydrodynamic state (flywheel pitch), we cannot use the `PTO` module we have used in previous tutorials. Instead, we will make two adjustments to our model compared to previous tutorials:\n", + "The generator shaft velocity is the difference between the buoy velocity and the flywheel velocity.\n", + "Because the velocity PTO depends not only on the hydrodynamic states (`x_wec`, buoy pitch) but also on an additional non-hydrodynamic state (flywheel pitch), we cannot use the `PTO` module we have used in previous tutorials.\n", + "We will make two adjustments to our previous modeling approach model to accurately model the PTO:\n", "\n", - "1. We will adjust the `x_opt` vector to include the flywheel motion\n", - "2. We will manually define the equations for the PTO dynamics and energy equations" + "1. Adjust the `x_opt` vector to include the flywheel motion\n", + "2. Manually define the PTO dynamics and energy equations in place on the `PTO` module" ] }, { @@ -306,11 +322,17 @@ "metadata": {}, "source": [ "#### Capturing flywheel motions\n", - "In previous tutorials, all the relevant WEC dynamics have been captured in the `x_wec` vector. However, `x_wec` assumes all the degrees of freedom are excited by waves, which is not the case for the flywheel in the Pioneer model.\n", + "In previous tutorials, all the relevant WEC dynamics have been captured in the `x_wec` vector.\n", + "However, `x_wec` assumes all the degrees of freedom are excited by waves, which is not the case for the flywheel in the Pioneer model.\n", + "Instead, we will include the flywheel dynamics in the `x_opt` vector.\n", "\n", - "The `x_opt` vector has previously only been used to model the PTO force. However, there is no limit to what `x_opt` can include; it can theoretically include anything that is needed to model the PTO control. Since we need to capture the flywheel dynamics to properly model the PTO, we will append the Fourier coefficients of the flywheel's position to `x_opt`. The first part of `x_opt` will continue to be used for the PTO force (e.g. Fourier coefficients of the force time-series for an unstructured controller). \n", + "The `x_opt` vector has previously only been used to model the PTO force.\n", + "However, there is no limit to what `x_opt` can include; it can theoretically include anything that is needed to model the PTO control.\n", + "Here, we will append the Fourier coefficients of the flywheel's position to `x_opt`.\n", + "The first part of `x_opt` will continue to be used for the PTO force (i.e. Fourier coefficients of the force time-series for an unstructured controller). \n", "\n", - "Note in the code further down that `x_opt` will be spliced as either `x_opt[:nstate_pto]` to exclude the flywheel position states, or `x_opt[nstate_pto:]` to exclude the PTO force states. Remember, `x_opt` can be defined arbitrarily, and specific meanings to the states are given elsewhere in the model." + "Note in the code further down that `x_opt` will be spliced as either `x_opt[:nstate_pto]` to exclude the flywheel position states, or `x_opt[nstate_pto:]` to exclude the PTO force states.\n", + "Remember, `x_opt` can be defined arbitrarily, and specific meanings to the states are given elsewhere in the model." ] }, { @@ -329,7 +351,8 @@ "metadata": {}, "source": [ "#### Manually defining PTO equations\n", - "Instead of calling the `PTO` module, we will manually define the dynamics and energy equations needed for WecOptTool to calculate electrical power, our quantity of interest for this optimization problem. These equations are also needed for the additional forces and constraints defined later on.\n", + "Instead of calling the `PTO` module, we will manually define the dynamics and energy equations needed for WecOptTool to calculate electrical power, our quantity of interest for this optimization problem.\n", + "These equations are also needed for the additional forces and constraints defined later on.\n", "\n", "##### Relative motion\n", "Here we define functions for the relative motion of the buoy and flywheel, derived from `x_wec` and `x_opt` respectively." @@ -383,7 +406,8 @@ "metadata": {}, "source": [ "##### PTO Impedance\n", - "The PTO impedance is defined using the same 2-port impedance model as in previous tutorials. The drivetrain inertia, friction, and stiffness are not included here since they are accounted for in the modeling of the additional non-hydrodynamic degree of freedom, as additional forces and constraints below." + "The PTO impedance is defined using the same 2-port impedance model as in previous tutorials.\n", + "The drivetrain inertia, friction, and stiffness are not included here since they are accounted for in the modeling of the additional non-hydrodynamic degree of freedom via additional forces and constraints below." ] }, { @@ -406,7 +430,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "When using the `PTO` module, the 2x2 impedance matrix is passed to the `wot.pto.PTO` object and converted into a transfer matrix to calculate the power at each frequency. While we're not using the `PTO` module, we will repeat this process almost verbatim here; see the `_make_abcd` and `_make_mimo_transfer_mat` functions in the `PTO` module if you are interested in more technical details." + "When using the `PTO` module, the 2x2 impedance matrix is passed to the `wot.pto.PTO` object and converted into a transfer matrix to calculate the power at each frequency.\n", + "While we are not using the `PTO` module, we will repeat this process almost verbatim here.\n", + "See the `_make_abcd` and `_make_mimo_transfer_mat` functions in the `PTO` module if you are interested in more technical details." ] }, { @@ -448,7 +474,8 @@ "metadata": {}, "source": [ "##### Power and Energy\n", - "Finally, we will wrap the relative motions, force, and impedance functions together to calculate average electrical power. This is again very close to what is found in the `PTO` module, except `rel_velocity` replaces `wot.pto.PTO.velocity` and `f_motor` replaces `wot.pto.PTO.force`." + "Finally, we will wrap the relative motions, force, and impedance functions together to calculate average electrical power.\n", + "This is again very close to what is found in the `PTO` module, except `rel_velocity` replaces `wot.pto.PTO.velocity` and `f_motor` replaces `wot.pto.PTO.force`." ] }, { @@ -493,7 +520,9 @@ "source": [ "### 1.4 Constraints\n", "The Pioneer model only contains one constraint based on the motor being used:\n", - " * **Peak torque** - The motor should avoid torques greater than 120 N-m. This is the same basic constraint as `const_f_pto` from Tutorial 1, except we use our `f_motor` definition instead of `wot.pto.PTO.force_on_wec`." + "\n", + " * **Peak torque** - The motor should avoid torques greater than 120 N-m.\n", + " This is the same basic constraint as `const_f_pto` from Tutorial 1, except we use our `f_motor` definition instead of `wot.pto.PTO.force_on_wec`." ] }, { @@ -514,9 +543,15 @@ "source": [ "### 1.5 Additional forces\n", "Here we add in the forces acting on the bodies other than those generated from wave loading.\n", - " * **Motor damping** - The resisting torque resulting from the damping on the PTO. This is very similar to the `wot.pto.PTO.force_on_wec` additional force given in previous tutorials. This is already being accounted for by the `f_motor` function above.\n", - " * **Flywheel friction** - The dissipative torque from the generator and components. This is defined as a nonlinear force using both Coulomb friction (dependent on the PTO direction) and the viscous friction (dependent on the PTO speed).\n", - " * **Magnetic spring** - The restoring torque caused by the torsional spring between the buoy and flywheel. Note the gear ratio is included in this equation, as the gearbox is connected in series with the spring. The position of the spring is the relative position of the flywheel and buoy scaled by the reciprocal of gear ratio squared. " + "\n", + " * **Motor damping** - The resisting torque resulting from the damping on the PTO.\n", + " This is very similar to the `wot.pto.PTO.force_on_wec` additional force given in previous tutorials. \n", + " This is already being accounted for by the `f_motor` function above.\n", + " * **Flywheel friction** - The dissipative torque from the generator and components.\n", + " This is defined as a nonlinear force using both Coulomb friction (dependent on the PTO direction) and the viscous friction (dependent on the PTO speed).\n", + " * **Magnetic spring** - The restoring torque caused by the torsional spring between the buoy and flywheel.\n", + " Note the gear ratio is included in this equation, as the gearbox is connected in series with the spring.\n", + " The position of the spring is the relative position of the flywheel and buoy scaled by the reciprocal of gear ratio squared. " ] }, { @@ -543,9 +578,15 @@ "metadata": {}, "source": [ "### 1.6 Flywheel residual equation\n", - "To make sure the flywheel dynamics are properly modeled in the pseudo-spectral method, we will introduce an equality constraint consisting of Newton's 2nd law of motion on the flywheel in residual form (i.e. $r(x) = I \\alpha - \\tau = 0 $). This is the same structure as the residual for the WEC dynamics described in the [WecOptTool theory documentation](https://sandialabs.github.io/WecOptTool/theory.html); also compare this equation to the 2nd dynamics equation listed [above](#pitch-resonator).\n", + "To make sure the flywheel dynamics are properly modeled in the pseudo-spectral method, we will introduce an equality constraint consisting of Newton's 2nd law of motion on the flywheel in residual form:\n", + "\n", + "$$ r(x) = I \\alpha - \\tau = 0 $$\n", "\n", - "This is required because the equations of motion described by `x_wec` only capture the buoy pitch degree of freedom. This provides a second equality constraint which will result in the correct two coupled dynamic equations." + "This is the same structure as the residual for the WEC dynamics described in the [WecOptTool theory documentation](https://sandialabs.github.io/WecOptTool/theory.html).\n", + "Also compare this equation to the 2nd dynamics equation listed [above](#pitch-resonator).\n", + "\n", + "This is required because the equations of motion described by `x_wec` only capture the buoy pitch degree of freedom.\n", + "This equation provides a second equality constraint to capture the flywheel pitch degree of freedom, resulting in two coupled dynamic equations." ] }, { @@ -573,7 +614,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now create the additional forces and constaints that will be passed to the `WEC` object. Note that `f_add` only pertains to the buoy since that is the only degree of freedom included in the BEM data; our flywheel residual equation handles these forces for the flywheel in the opposite direction." + "We now create the additional forces and constaints that will be passed to the `WEC` object.\n", + "Note that `f_add` only pertains to the buoy since that is the only degree of freedom included in the BEM data; our flywheel residual equation handles these forces for the flywheel in the opposite direction." ] }, { @@ -620,8 +662,9 @@ "metadata": {}, "source": [ "## 2. Regular wave results\n", + "\n", "### 2.1 Solve\n", - "As always, we will optimize for electrical power absorption." + "As in previous tutorials, we will optimize for electrical power absorption." ] }, { @@ -647,7 +690,9 @@ "metadata": {}, "source": [ "### 2.2 Post-process and plot\n", - "Again, since we're not using the `PTO` module, post-processing using `wot.pto.PTO.post_process` is not an option here, so we have to manually post-process the outputs related to the PTO and flywheel. This is pretty intuitive using the functions we created earlier. The outputs related to the buoy can still be derived directly from `wot.wec.post_process`." + "Again, since we are not using the `PTO` module, post-processing using `wot.pto.PTO.post_process` is not an option here, so we have to manually post-process the outputs related to the PTO and flywheel.\n", + "This is pretty intuitive using the functions we created earlier.\n", + "The outputs related to the buoy can still be derived directly from `wot.wec.post_process`." ] }, { @@ -657,7 +702,7 @@ "outputs": [], "source": [ "nsubsteps = 5\n", - "wec_fdom, wec_tdom = wec.post_process(results[0], waves_regular.sel(realization=0), nsubsteps=nsubsteps)\n", + "wec_fdom, wec_tdom = wec.post_process(wec, results, waves_regular, nsubsteps=nsubsteps)\n", "\n", "# Manually post-process PTO and flywheel outputs\n", "x_wec, x_opt = wot.decompose_state(results[0].x, 1, nfreq)\n", @@ -679,7 +724,9 @@ "Some observations about the optimized Pioneer model:\n", "\n", " * The flywheel and buoy match frequency, but are out of phase due to the forcing from the spring and motor applied to both bodies.\n", - " * The gearing in series with the torsional spring (which reduces the effective spring stiffness) significantly amplifies the rotation of the flywheel despite the buoy only pitching modestly. The gear ratio of 3 was selected by the Pioneer team to keep the system in resonance considering the spring stiffness, moment of inertia of the flywheel, and the resonance frequency of the buoy. Try increasing the gear ratio or flywheel moment of inertia -- this will actually reduce the power generated since it will bring the effective stiffness of the flywheel out of resonance, causing less motion relative to the buoy.\n", + " * The gearing in series with the torsional spring (which reduces the effective spring stiffness) significantly amplifies the rotation of the flywheel despite the buoy only pitching modestly.\n", + " The gear ratio of 3 was selected by the Pioneer team to keep the system in resonance considering the spring stiffness, moment of inertia of the flywheel, and the resonance frequency of the buoy.\n", + " Try increasing the gear ratio or flywheel moment of inertia—this will actually reduce the power generated since it will bring the effective stiffness of the flywheel out of resonance, causing less motion relative to the buoy.\n", " * The mechanical and electrical power outputs are similar, since our impedance model has only a small amount of resistance and no inductance.\n", " * The buoy's pitch amplitude is larger than expected for this device and is likely due to underestimation of the radiation damping by the BEM. " ] @@ -691,7 +738,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=4, sharex=True, figsize=(12, 12))\n", - "t = wec_tdom.time.values\n", + "t = wec_tdom[0].time.values\n", "\n", "# Positions\n", "ax[0].plot(t, fw_pos*180/np.pi, label='Flywheel', c='k')\n", @@ -702,7 +749,7 @@ "\n", "ax0 = ax[0].twinx()\n", "ax0.tick_params(axis='y', labelcolor='b')\n", - "(wec_tdom.pos*180/np.pi).plot(hue='influenced_dof', label='Buoy', ax=ax0, c='b')\n", + "(wec_tdom[0].pos*180/np.pi).plot(hue='influenced_dof', label='Buoy', ax=ax0, c='b')\n", "ax0.set_ylabel('Buoy pos. [deg]', color='b')\n", "ax0.set_title('')\n", "\n", @@ -715,7 +762,7 @@ "\n", "ax1 = ax[1].twinx()\n", "ax1.tick_params(axis='y', labelcolor='b')\n", - "wec_tdom['force'].sel(type=['Froude_Krylov', 'diffraction']).sum('type').plot(ax=ax1, c='b')\n", + "wec_tdom[0]['force'].sel(type=['Froude_Krylov', 'diffraction']).sum('type').plot(ax=ax1, c='b')\n", "ax1.set_ylabel('Excitation torque [N-m]', color='b')\n", "ax1.set_title('Torque')\n", "\n", @@ -742,14 +789,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 2.3 Sankey Diagram" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We use a Sankey diagram to visualize the average power flow from waves to electricity. On the very left is the theoretically possible optimal excitation. In order to reach this upper bound of exitation the WEC pitch velocity would need to be in phase with the pitch excitation force, in this case the radiated power is equal the amount of absorbed mechanical power. In practice, this will usually imply putting electrical power into the system (something we want to avoid!!!). \n", + "### 2.3 Sankey Diagram\n", + "We use a Sankey diagram to visualize the average power flow from waves to electricity.\n", + "The left end shows the theoretically possible optimal excitation.\n", + "In order to reach this upper bound of excitation, the WEC pitch velocity would need to be in phase with the pitch excitation force (or, in this case, the radiated power is equal the amount of absorbed mechanical power).\n", + "In practice, this will usually imply putting electrical power into the system (something we want to avoid!!!).\n", + "\n", "With co-design, we are instead tapping the unused potential while limiting PTO losses and radiated power." ] }, @@ -760,26 +805,34 @@ "outputs": [], "source": [ "hydro_data = wot.add_linear_friction(bem_data, friction=None) \n", - "hydro_data = wot.check_linear_damping(hydro_data, uniform_shift=False) \n", + "hydro_data = wot.check_radiation_damping(hydro_data, uniform_shift=False) \n", "\n", "Zi = wot.hydrodynamic_impedance(hydro_data)\n", "Rad_res = np.real(Zi.squeeze())\n", "\n", - "Fex = wec_fdom.force.sel(type=['Froude_Krylov', 'diffraction']).sum('type')\n", - "Vel = wec_fdom.vel\n", + "Fex = wec_fdom[0].force.sel(type=['Froude_Krylov', 'diffraction']).sum('type')\n", + "Vel = wec_fdom[0].vel\n", "\n", "P_max_absorbable = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item() # after Falnes Eq. 6.10\n", "P_opt_excitation = 2*P_max_absorbable # after Falnes Eq. 6.10\n", "P_radiated = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item()) # after Falnes Eq. 6.4\n", "P_excited= (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item() # after Falnes Eq. 6.3\n", - "P_absorbed = P_excited - P_radiated # after Falnes Eq. 6.2 absorbed by WEC-PTO (difference between actual excitation power and actual radiated power needs to be absorbed by PTO)" + "\n", + "power_flows = {'Optimal Excitation' : -2* (np.real(P_max_absorbable)), \n", + " 'Radiated': -1*(np.real(P_radiated)), \n", + " 'Actual Excitation': -1*(np.real(P_excited)), \n", + "}\n", + "power_flows['Unused Potential'] = power_flows['Optimal Excitation'] - power_flows['Actual Excitation']\n", + "power_flows['Absorbed'] = power_flows['Actual Excitation'] - power_flows['Radiated'] # after Falnes Eq. 6.2 \n", + "#absorbed by WEC-PTO (difference between actual excitation power and actual radiated power needs to be absorbed by PTO)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We also calculate the power dissipated due to flywheel friction and make sure that the absorbed power (calculated as the difference between excited and radiated power) matches the sum of mechanical power captured by the PTO and the power dissipated due to flywheel friction. We use a relative tolerance of 1%." + "We also calculate the power dissipated due to flywheel friction and make sure that the absorbed power (calculated as the difference between excited and radiated power) matches the sum of mechanical power captured by the PTO and the power dissipated due to flywheel friction.\n", + "We use a relative tolerance of 1%." ] }, { @@ -800,9 +853,7 @@ " return vel_td * force_td\n", "\n", "fw_fric_power = fw_friction_power(wec, x_wec, x_opt, waves_regular, nsubsteps)\n", - "avg_fw_fric_power = np.mean(fw_fric_power)\n", - "\n", - "assert(np.isclose(P_absorbed, -1*(avg_mech_power -avg_fw_fric_power), rtol=0.01)) # assert that solver solution matches" + "avg_fw_fric_power = np.mean(fw_fric_power)\n" ] }, { @@ -811,39 +862,18 @@ "metadata": {}, "outputs": [], "source": [ - "from matplotlib.sankey import Sankey\n", - "\n", - "P_PTO_loss = avg_mech_power - avg_elec_power \n", - "P_unused = P_opt_excitation - P_excited # Difference between the theoretical optimum excitation, if the WEC velocity would be in resonance with the excitation force\n", - "\n", - "Power_flows = [P_opt_excitation, P_PTO_loss, -1*avg_fw_fric_power, -1*P_radiated, -1*P_unused, avg_elec_power, ]\n", - "\n", - "fig = plt.figure(figsize = [6,4])\n", - "ax = fig.add_subplot(1, 1, 1,)\n", - "sankey = Sankey(ax=ax, \n", - " scale=0.5/P_max_absorbable,\n", - " offset= 0,\n", - " format = '%.2f W',shoulder = 0.02)\n", - "\n", - "sankey.add(flows=Power_flows, \n", - " labels = ['Optimal Excitation \\n $ 2 \\\\frac{\\left | F_e \\\\right |^2}{8Re(Z_i)} = 2*P_{mech}^{max}$ ', \n", - " 'PTO-Loss \\n $ P_{mech} - P_{elec}$', \n", - " 'Flywheel friction',\n", - " 'Radiated \\n $ \\\\frac{1}{2} Re(Z_i) \\left | U \\\\right |^2 $ ', \n", - " 'Unused Potential \\n(neither absorbed nor radiated)', \n", - " 'Electrical'], \n", - " orientations=[0, -1, -1,-1, -1, 0], # arrow directions\n", - " pathlengths = [0.4,0.2,0.6,0.6,0.7,0.5,],\n", - " trunklength = 1.5,\n", - " edgecolor = '#2b8cbe',\n", - " facecolor = '#2b8cbe' )\n", - "\n", - "diagrams = sankey.finish()\n", - "for diagram in diagrams:\n", - " for text in diagram.texts:\n", - " text.set_fontsize(10);\n", - "plt.axis(\"off\") \n", - "plt.show()" + "power_flows['Electrical (solver)']= avg_elec_power #solver determines sign\n", + "power_flows['Mechanical (solver)']= avg_mech_power - avg_fw_fric_power #solver determines sign, friction is dissipated\n", + "power_flows['PTO Loss'] = power_flows['Mechanical (solver)'] - power_flows['Electrical (solver)'] " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "wot.utilities.plot_power_flow(power_flows)" ] }, { @@ -851,10 +881,16 @@ "metadata": {}, "source": [ "## 3. Irregular wave\n", + "\n", "### 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. " + "Interestingly, due to the narrow banded resonance of the flywheel, the controller attempts to make the excitation force monochromatic with the resonant frequency.\n", + "To achieve this, it uses significant reactive power (power by the PTO into the system).\n", + "This is still worth it though, resulting in a larger average electrical power output.\n", + "\n", + "As noted previously, the optimization problem is solved for each wave realization.\n", + "The optimal average power shown is the total average across the different realizations." ] }, { @@ -890,7 +926,7 @@ "outputs": [], "source": [ "nsubsteps = 5\n", - "wec_fdom, wec_tdom = wec.post_process(results[0], waves_irregular.sel(realization=0), nsubsteps=nsubsteps)\n", + "wec_fdom, wec_tdom = wec.post_process(wec,results, waves_irregular, nsubsteps=nsubsteps)\n", "\n", "# Manually post-process PTO and flywheel outputs\n", "x_wec, x_opt = wot.decompose_state(results[0].x, 1, nfreq)\n", @@ -912,7 +948,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=4, sharex=True, figsize=(12, 12))\n", - "t = wec_tdom.time.values\n", + "t = wec_tdom[0].time.values\n", "\n", "# Positions\n", "ax[0].plot(t, fw_pos*180/np.pi, label='Flywheel', c='k')\n", @@ -923,7 +959,7 @@ "\n", "ax0 = ax[0].twinx()\n", "ax0.tick_params(axis='y', labelcolor='b')\n", - "(wec_tdom.pos*180/np.pi).plot(hue='influenced_dof', label='Buoy', ax=ax0, c='b')\n", + "(wec_tdom[0].pos*180/np.pi).plot(hue='influenced_dof', label='Buoy', ax=ax0, c='b')\n", "ax0.set_ylabel('Buoy pos. [deg]', color='b')\n", "ax0.set_title('')\n", "\n", @@ -936,7 +972,7 @@ "\n", "ax1 = ax[1].twinx()\n", "ax1.tick_params(axis='y', labelcolor='b')\n", - "wec_tdom['force'].sel(type=['Froude_Krylov', 'diffraction']).sum('type').plot(ax=ax1, c='b')\n", + "wec_tdom[0]['force'].sel(type=['Froude_Krylov', 'diffraction']).sum('type').plot(ax=ax1, c='b')\n", "ax1.set_ylabel('Excitation torque [N-m]', color='b')\n", "ax1.set_title('Torque')\n", "\n", @@ -959,12 +995,131 @@ " axi.autoscale(axis='x', tight=True)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.3 Notes on selection of frequency array and realizations\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.\n", + "Often, capturing nonlinear dynamics (such as due to constraints) requires a larger frequency range.\n", + "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).\n", + "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.\n", + "The frequency values are normalized in the plots below to more clearly show the harmonics.\n", + "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.\n", + "The Pioneer WEC has a very narrow banded response, which means it requires a large number of frequencies to model the response accurately.\n", + "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.\n", + "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).\n", + "As shown below, an array of 150 frequencies is sufficient to model the system to within about 2% of the actual resultant mean power.\n", + "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": [] + "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.\n", + "Generally, it is recommended that enough realizations be used for the total simulation time to equal 20 minutes. \n", + "\n", + "$$ t_{total} = \\frac{n_{realizations}}{f1} $$\n", + "\n", + "Although it usually leads to consistent results, this recommendation is somewhat arbitrary.\n", + "To better understand the effect of the number of realizations on the overall result, it can be useful to complete a convergence study.\n", + "The actual number of realizations needed depends on the problem itself (dynamics, constraints, etc.) and the desired precision of the result.\n", + "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]')" + ] } ], "metadata": { @@ -983,7 +1138,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.11.4" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index af0cc9737..34fdb8da7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "wecopttool" -version = "2.7.0" +version = "3.0.0" description = "WEC Design Optimization Toolbox" readme = "README.md" authors = [ @@ -20,13 +20,13 @@ classifiers = [ ] dependencies = [ "numpy>=1.20", - "scipy", "xarray", "autograd", "capytaine", "joblib", "wavespectra>=3.13", "netcdf4", + "cyipopt", ] [project.optional-dependencies] diff --git a/tests/test_core.py b/tests/test_core.py index b8beb4fee..55f093f0e 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -919,8 +919,8 @@ def test_read_write_netcdf(self, hydro_data, hydro_data_new): assert hydro_data.equals(hydro_data_new) -class TestCheckLinearDamping: - """Test function :python:`check_linear_damping`.""" +class TestCheckRadiationDamping: + """Test function :python:`check_radiation_damping`.""" @pytest.fixture(scope="class") def data(self, hydro_data): @@ -933,44 +933,39 @@ def data(self, hydro_data): @pytest.fixture(scope="class") def tol(self, data): - """Tolerance for function :python:`check_linear_damping`.""" + """Tolerance for function :python:`check_radiation_damping`.""" return 0.01 @pytest.fixture(scope="class") def data_new_uniform(self, data, tol): """Hydrodynamic data structure for which the function - :python:`check_linear_damping` has been called. + :python:`check_radiation_damping` has been called. """ - return wot.check_linear_damping(data, tol) + return wot.check_radiation_damping(data, tol, True) @pytest.fixture(scope="class") def data_new_nonuniform(self, data, tol): """Hydrodynamic data structure for which the function - :python:`check_linear_damping` has been called. + :python:`check_radiation_damping` has been called. """ - return wot.check_linear_damping(data, tol, False) + return wot.check_radiation_damping(data, tol, False) - def test_friction(self, data_new_uniform, tol): + def test_radiation(self, data, data_new_uniform, tol): """Test that the modified friction diagonal has the expected value for a uniform shift. """ - assert np.allclose(np.diagonal(data_new_uniform.friction.values), tol) + radiation_diff = (data_new_uniform.radiation_damping.values + - data.radiation_damping.values) + assert np.allclose(radiation_diff, radiation_diff[0], tol) - def test_only_diagonal_friction(self, data, data_new_uniform): - """Test that only the diagonal was changed for a uniform shift.""" - data_org = data.copy(deep=True) - def nodiag(x): - return x.friction.values - np.diag(np.diagonal(x.friction.values)) - assert np.allclose(nodiag(data_new_uniform), nodiag(data_org)) - - def test_only_friction(self, data, data_new_uniform): - """Test that only the friction is changed in the hydrodynamic + def test_only_radiation(self, data, data_new_uniform): + """Test that only the radiation is changed in the hydrodynamic data for a uniform shift. """ - data_new_nofric = data_new_uniform.copy(deep=True - ).drop_vars('friction') - data_org_nofric = data.copy(deep=True).drop_vars('friction') - assert data_new_nofric.equals(data_org_nofric) + data_new_norad = data_new_uniform.copy(deep=True + ).drop_vars('radiation_damping') + data_org_norad = data.copy(deep=True).drop_vars('radiation_damping') + assert data_new_norad.equals(data_org_norad) def test_damping(self, data_new_nonuniform, tol): """Test that the modified radiation damping diagonal has the expected diff --git a/tests/test_integration.py b/tests/test_integration.py index ea116fa40..0e9bc6dea 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -10,7 +10,7 @@ kplim = -1e1 - +min_damping = 45 @pytest.fixture(scope="module") def f1(): @@ -73,10 +73,13 @@ def fb(): @pytest.fixture(scope="module") -def bem(f1, nfreq, fb): - """Boundary elemement model (Capytaine) results""" +def hydro_data(f1, nfreq, fb): + """Boundary element model (Capytaine) results (with friction added)""" freq = wot.frequency(f1, nfreq, False) - return wot.run_bem(fb, freq) + hydro_data = wot.run_bem(fb, freq) + hd = wot.add_linear_friction(hydro_data) + hd = wot.check_radiation_damping(hd, min_damping=min_damping) + return hd @pytest.fixture(scope='module') @@ -91,21 +94,26 @@ def regular_wave(f1, nfreq): @pytest.fixture(scope='module') -def irregular_wave(f1, nfreq): +def long_crested_wave(f1, nfreq): """Idealized (Pierson-Moskowitz) spectrum wave""" freq = wot.frequency(f1, nfreq, False) fp = 0.3 hs = 0.0625*1.9 - spec = wot.waves.pierson_moskowitz_spectrum(freq, fp, hs) - waves = wot.waves.long_crested_wave(spec) + spec_fun = lambda f: wot.waves.pierson_moskowitz_spectrum(freq=f, + fp=fp, + hs=hs) + efth = wot.waves.omnidirectional_spectrum(f1=f1, nfreq=nfreq, + spectrum_func=spec_fun, + ) + waves = wot.waves.long_crested_wave(efth, nrealizations=2) return waves @pytest.fixture(scope='module') -def wec_from_bem(f1, nfreq, bem, fb, pto): +def wec_from_bem(f1, nfreq, hydro_data, fb, pto): """Simple WEC: 1 DOF, no constraints.""" f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_bem(bem, f_add=f_add) + wec = wot.WEC.from_bem(hydro_data, f_add=f_add) return wec @@ -113,14 +121,15 @@ def wec_from_bem(f1, nfreq, bem, fb, pto): def wec_from_floatingbody(f1, nfreq, fb, pto): """Simple WEC: 1 DOF, no constraints.""" f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add) + wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add, + min_damping=min_damping) return wec @pytest.fixture(scope='module') -def wec_from_impedance(bem, pto, fb): +def wec_from_impedance(hydro_data, pto, fb): """Simple WEC: 1 DOF, no constraints.""" - bemc = bem.copy() + bemc = hydro_data.copy() omega = bemc['omega'].values w = np.expand_dims(omega, [1,2]) A = bemc['added_mass'].values @@ -134,45 +143,24 @@ def wec_from_impedance(bem, pto, fb): freqs = omega / (2 * np.pi) impedance = (A + mass)*(1j*w) + B + K/(1j*w) - exc_coeff = bem['Froude_Krylov_force'] + bem['diffraction_force'] + exc_coeff = hydro_data['Froude_Krylov_force'] + hydro_data['diffraction_force'] f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add) + wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add, + min_damping=min_damping) return wec @pytest.fixture(scope='module') -def resonant_wave(f1, nfreq, fb, bem): +def resonant_wave(f1, nfreq, fb, hydro_data): """Regular wave at natural frequency of the WEC""" - hd = wot.add_linear_friction(bem) + hd = wot.add_linear_friction(hydro_data) Zi = wot.hydrodynamic_impedance(hd) wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item() waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1) return waves -def test_solve_callback(wec_from_bem, regular_wave, pto, nfreq, capfd): - """Check that user can set a custom callback""" - - cbstring = 'hello world!' - - def my_callback(my_wec, x_wec, x_opt, wave): - print(cbstring) - - _ = wec_from_bem.solve(regular_wave, - obj_fun=pto.average_power, - nstate_opt=2*nfreq, - scale_x_wec=1.0, - scale_x_opt=0.01, - scale_obj=1e-1, - callback=my_callback, - optim_options={'maxiter': 1}) - - out, err = capfd.readouterr() - - assert out.split('\n')[0] == cbstring - - @pytest.mark.parametrize("bounds_opt", [Bounds(lb=kplim, ub=0), ((kplim, 0),)]) def test_solve_bounds(bounds_opt, wec_from_bem, regular_wave, @@ -231,15 +219,48 @@ def hstiff(self, fb): return fb.compute_hydrostatic_stiffness() @pytest.fixture(scope='class') - def hydro_impedance(self, bem): + def hydro_impedance(self, hydro_data): """Intrinsic hydrodynamic impedance""" - hd = wot.add_linear_friction(bem) - hd = wot.check_linear_damping(hd) - Zi = wot.hydrodynamic_impedance(hd) + Zi = wot.hydrodynamic_impedance(hydro_data) return Zi + + @pytest.fixture(scope='class') + def unstruct_wec(self, + hydro_data, + pto): + """WaveBot WEC object with unstructured controller""" + + f_add = {"PTO": pto.force_on_wec} + wec = wot.WEC.from_bem(hydro_data, f_add=f_add) + + return wec + + @pytest.fixture(scope='class') + def long_crested_wave_unstruct_res(self, + unstruct_wec, + long_crested_wave, + pto, + hydro_data, + nfreq): + """Solution for an unstructured controller with multiple long crested + waves""" + + f_add = {"PTO": pto.force_on_wec} + wec = wot.WEC.from_bem(hydro_data, f_add=f_add) + + res = unstruct_wec.solve(waves=long_crested_wave, + obj_fun=pto.average_power, + nstate_opt=2*nfreq, + x_wec_0=1e-3*np.ones(wec.nstate_wec), + scale_x_wec=1e1, + scale_x_opt=1e-3, + scale_obj=5e-2, + ) + + return res def test_p_controller_resonant_wave(self, - bem, + hydro_data, resonant_wave, p_controller_pto, hydro_impedance): @@ -247,7 +268,7 @@ def test_p_controller_resonant_wave(self, wave""" f_add = {"PTO": p_controller_pto.force_on_wec} - wec = wot.WEC.from_bem(bem, f_add=f_add) + wec = wot.WEC.from_bem(hydro_data, f_add=f_add) res = wec.solve(waves=resonant_wave, obj_fun=p_controller_pto.average_power, @@ -257,14 +278,14 @@ def test_p_controller_resonant_wave(self, scale_x_wec=1e2, scale_x_opt=1e-3, scale_obj=1e-1, - optim_options={'ftol': 1e-10}, + optim_options={'tol': 1e-10}, bounds_opt=((-1*np.infty, 0),), ) power_sol = -1*res[0]['fun'] - res_fd, _ = wec.post_process(res[0], resonant_wave.sel(realization=0), nsubsteps=1) - Fex = res_fd.force.sel( + res_fd, _ = wec.post_process(wec, res, resonant_wave, nsubsteps=1) + Fex = res_fd[0].force.sel( type=['Froude_Krylov', 'diffraction']).sum('type') power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze()) ).squeeze().sum('omega').item() @@ -272,14 +293,14 @@ def test_p_controller_resonant_wave(self, assert power_sol == approx(power_optimal, rel=0.03) def test_pi_controller_regular_wave(self, - bem, + hydro_data, regular_wave, pi_controller_pto, hydro_impedance): """PI controller matches optimal for any regular wave""" f_add = {"PTO": pi_controller_pto.force_on_wec} - wec = wot.WEC.from_bem(bem, f_add=f_add) + wec = wot.WEC.from_bem(hydro_data, f_add=f_add) res = wec.solve(waves=regular_wave, obj_fun=pi_controller_pto.average_power, @@ -295,47 +316,46 @@ def test_pi_controller_regular_wave(self, power_sol = -1*res[0]['fun'] - res_fd, _ = wec.post_process(res[0], regular_wave.sel(realization=0), nsubsteps=1) - Fex = res_fd.force.sel( + res_fd, _ = wec.post_process(wec, res, regular_wave, nsubsteps=1) + Fex = res_fd[0].force.sel( type=['Froude_Krylov', 'diffraction']).sum('type') power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze()) ).squeeze().sum('omega').item() assert power_sol == approx(power_optimal, rel=1e-4) - def test_unstructured_controller_irregular_wave(self, - fb, - bem, - regular_wave, - pto, - nfreq, - hydro_impedance): + + def test_unstructured_controller_long_crested_wave(self, + unstruct_wec, + long_crested_wave, + hydro_impedance, + long_crested_wave_unstruct_res, + pto): """Unstructured (numerical optimal) controller matches optimal for any - irregular wave when unconstrained""" - - f_add = {"PTO": pto.force_on_wec} - wec = wot.WEC.from_bem(bem, f_add=f_add) - - res = wec.solve(waves=regular_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, - ) + irregular (long crested) wave when unconstrained""" - power_sol = -1*res[0]['fun'] + power_sol = -1*long_crested_wave_unstruct_res[0]['fun'] - res_fd, _ = wec.post_process(res[0], regular_wave.sel(realization=0), nsubsteps=1) - Fex = res_fd.force.sel( + res_fd, _ = unstruct_wec.post_process(unstruct_wec, long_crested_wave_unstruct_res, + long_crested_wave, + nsubsteps=1) + Fex = res_fd[0].force.sel( type=['Froude_Krylov', 'diffraction']).sum('type') power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze()) ).squeeze().sum('omega').item() - assert power_sol == approx(power_optimal, rel=1e-3) + assert power_sol == approx(power_optimal, rel=1e-2) + + def test_unconstrained_solutions_multiple_phase_realizations(self, + long_crested_wave_unstruct_res): + """Solutions for average power with an unstructured controller + (no constraints) match for different phase realizations""" + + pow = [res['fun'] for res in long_crested_wave_unstruct_res] + + assert pow[0] == approx(pow[1], rel=1e-4) def test_saturated_pi_controller(self, - bem, + hydro_data, regular_wave, pto, nfreq): @@ -357,7 +377,7 @@ def const_f_pto(wec, x_wec, x_opt, waves): f = pto['us'].force_on_wec(wec, x_wec, x_opt, waves, nsubsteps=4) return f_max - np.abs(f.flatten()) - wec['us'] = wot.WEC.from_bem(bem, + wec['us'] = wot.WEC.from_bem(hydro_data, f_add={"PTO": pto['us'].force_on_wec}, constraints=[{'type': 'ineq', 'fun': const_f_pto, }]) @@ -372,7 +392,7 @@ def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1): pto['pi'] = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof), controller=saturated_pi,) - wec['pi'] = wot.WEC.from_bem(bem, + wec['pi'] = wot.WEC.from_bem(hydro_data, f_add={"PTO": pto['pi'].force_on_wec}, constraints=[]) @@ -405,13 +425,13 @@ def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1): nsubstep_postprocess = 4 pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key], - res[key][0], - regular_wave.sel(realization=0), + res[key], + regular_wave, nsubstep_postprocess) - xr.testing.assert_allclose(pto_tdom['pi'].power.squeeze().mean('time'), - pto_tdom['us'].power.squeeze().mean('time'), + xr.testing.assert_allclose(pto_tdom['pi'][0].power.squeeze().mean('time'), + pto_tdom['us'][0].power.squeeze().mean('time'), rtol=1e-1) - xr.testing.assert_allclose(pto_tdom['us'].force.max(), - pto_tdom['pi'].force.max()) + xr.testing.assert_allclose(pto_tdom['us'][0].force.max(), + pto_tdom['pi'][0].force.max()) diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 000000000..f2a9f6be4 --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,209 @@ +""" Unit tests for functions in the :python:`utilities.py` module. +""" + +import pytest +import numpy as np +import xarray as xr +from matplotlib.pyplot import Figure, Axes +import wecopttool as wot +from pytest import approx +import capytaine as cpy + + + +# test function in the utilities.py + + + +@pytest.fixture(scope="module") +def power_flows(): + """Dictionary of power flows.""" + pflows = {'Optimal Excitation': -100, + 'Radiated': -20, + 'Actual Excitation': -70, + 'Electrical (solver)': -40, + 'Mechanical (solver)': -50, + 'Absorbed': -50, + 'Unused Potential': -30, + 'PTO Loss': -10 + } + return pflows + +@pytest.fixture(scope="module") +def f1(): + """Fundamental frequency [Hz].""" + return 0.1 + + +@pytest.fixture(scope="module") +def nfreq(): + """Number of frequencies in frequency vector.""" + return 5 + +@pytest.fixture(scope="module") +def ndof(): + """Number of degrees of freedom.""" + return 2 + +@pytest.fixture(scope="module") +def ndir(): + """Number of wave directions.""" + return 3 + +@pytest.fixture(scope='module') +def bem_data(f1, nfreq, ndof, ndir): + """Synthetic BEM data.""" + # TODO - start using single BEM solution across entire test suite + coords = { + 'omega': [2*np.pi*(ifreq+1)*f1 for ifreq in range(nfreq)], + 'influenced_dof': [f'DOF_{idof+1}' for idof in range(ndof)], + 'radiating_dof': [f'DOF_{idof+1}' for idof in range(ndof)], + 'wave_direction': [2*np.pi/ndir*idir for idir in range(ndir)], + } + radiation_dims = ['omega', 'radiating_dof', 'influenced_dof'] + excitation_dims = ['omega', 'influenced_dof', 'wave_direction'] + hydrostatics_dims = ['radiating_dof', 'influenced_dof'] + + added_mass = np.ones([nfreq, ndof, ndof]) + radiation_damping = np.ones([nfreq, ndof, ndof]) + diffraction_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j + Froude_Krylov_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j + inertia_matrix = np.ones([ndof, ndof]) + hydrostatic_stiffness = np.ones([ndof, ndof]) + + data_vars = { + 'added_mass': (radiation_dims, added_mass), + 'radiation_damping': (radiation_dims, radiation_damping), + 'diffraction_force': (excitation_dims, diffraction_force), + 'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force), + 'inertia_matrix': (hydrostatics_dims, inertia_matrix), + 'hydrostatic_stiffness': (hydrostatics_dims, hydrostatic_stiffness) + } + return xr.Dataset(data_vars=data_vars, coords=coords) + +@pytest.fixture(scope='module') +def intrinsic_impedance(bem_data): + bem_data = wot.add_linear_friction(bem_data) + intrinsic_impedance = wot.hydrodynamic_impedance(bem_data) + return intrinsic_impedance + +@pytest.fixture(scope='module') +def pi_controller_pto(): + """Basic PTO: proportional-integral (PI) controller, 1DOF, mechanical + power.""" + ndof = 1 + pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof), + controller=wot.pto.controller_pi, + names=["PI controller PTO"]) + return pto + +@pytest.fixture(scope='module') +def regular_wave(f1, nfreq): + """Single frequency wave""" + wfreq = 0.3 + wamp = 0.0625 + wphase = 0 + wdir = 0 + waves = wot.waves.regular_wave(f1, nfreq, wfreq, wamp, wphase, wdir) + return waves + +@pytest.fixture(scope="module") +def fb(): + """Capytaine FloatingBody object""" + try: + import wecopttool.geom as geom + except ImportError: + pytest.skip( + 'Skipping integration tests due to missing optional geometry ' + + 'dependencies. Run `pip install wecopttool[geometry]` to run ' + + 'these tests.' + ) + mesh_size_factor = 0.5 + wb = geom.WaveBot() + mesh = wb.mesh(mesh_size_factor) + fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot") + fb.add_translation_dof(name="Heave") + return fb + + +@pytest.fixture(scope="module") +def wb_bem(f1, nfreq, fb): + """Boundary elemement model (Capytaine) results""" + freq = wot.frequency(f1, nfreq, False) + return wot.run_bem(fb, freq) + +@pytest.fixture(scope='class') +def wb_hydro_impedance(wb_bem): + """Intrinsic hydrodynamic impedance""" + hd = wot.add_linear_friction(wb_bem) + hd = wot.check_radiation_damping(hd) + Zi = wot.hydrodynamic_impedance(hd) + return Zi + + + + +def test_plot_hydrodynamic_coefficients(bem_data,ndof): + bem_figure_list = wot.utilities.plot_hydrodynamic_coefficients(bem_data) + correct_len = ndof*(ndof+1)/2 #using only the subdiagonal elements + #added mass + fig_am = bem_figure_list[0][0] + assert correct_len == len(fig_am.axes) + assert isinstance(fig_am,Figure) + #radiation damping + fig_rd = bem_figure_list[1][0] + assert correct_len == len(fig_rd.axes) + assert isinstance(fig_rd,Figure) + #radiation damping + fig_ex = bem_figure_list[2][0] + assert ndof == len(fig_ex.axes) + assert isinstance(fig_ex,Figure) + +def test_plot_bode_impedance(intrinsic_impedance, ndof): + fig_Zi, axes_Zi = wot.utilities.plot_bode_impedance(intrinsic_impedance) + + assert 2*ndof*ndof == len(fig_Zi.axes) + assert isinstance(fig_Zi,Figure) + assert all([isinstance(ax, Axes) for ax in np.reshape(axes_Zi,-1)]) + + +def test_plot_power_flow(power_flows): + fig_sankey, ax_sankey = wot.utilities.plot_power_flow(power_flows) + + assert isinstance(fig_sankey, Figure) + assert isinstance(ax_sankey, Axes) + +def test_calculate_power_flow(wb_bem, + regular_wave, + pi_controller_pto, + wb_hydro_impedance): + """PI controller matches optimal for any regular wave, + thus we check if the radiated power is equal the absorber power + and if the Optimal excitation is equal the actual excitation""" + + f_add = {"PTO": pi_controller_pto.force_on_wec} + wec = wot.WEC.from_bem(wb_bem, f_add=f_add) + + res = wec.solve(waves=regular_wave, + obj_fun=pi_controller_pto.average_power, + nstate_opt=2, + x_wec_0=1e-1*np.ones(wec.nstate_wec), + x_opt_0=[-1e3, 1e4], + scale_x_wec=1e2, + scale_x_opt=1e-3, + scale_obj=1e-2, + optim_options={'maxiter': 50}, + bounds_opt=((-1e4, 0), (0, 2e4),) + ) + + pflows = wot.utilities.calculate_power_flows(wec, + pi_controller_pto, + res, + regular_wave, + wb_hydro_impedance) + + assert pflows['Absorbed'] == approx(pflows['Radiated'], rel=1e-4) + assert pflows['Optimal Excitation'] == approx(pflows['Actual Excitation'], rel=1e-4) + + + diff --git a/wecopttool/__init__.py b/wecopttool/__init__.py index 9b3d6b666..3cd720358 100644 --- a/wecopttool/__init__.py +++ b/wecopttool/__init__.py @@ -32,6 +32,7 @@ from wecopttool.core import * from wecopttool import waves from wecopttool import pto +from wecopttool import utilities try: from wecopttool import geom diff --git a/wecopttool/core.py b/wecopttool/core.py index 96cf32a48..065327b91 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -31,7 +31,7 @@ "td_to_fd", "read_netcdf", "write_netcdf", - "check_linear_damping", + "check_radiation_damping", "check_impedance", "force_from_rao_transfer_function", "force_from_impedance", @@ -70,7 +70,7 @@ import capytaine as cpy from scipy.optimize import minimize, OptimizeResult, Bounds from scipy.linalg import block_diag, dft - +from cyipopt import minimize_ipopt # logger _log = logging.getLogger(__name__) @@ -289,7 +289,7 @@ def from_bem( f_add: Optional[TIForceDict] = None, constraints: Optional[Iterable[Mapping]] = None, min_damping: Optional[float] = _default_min_damping, - uniform_shift: Optional[bool] = True, + uniform_shift: Optional[bool] = False, dof_names: Optional[Iterable[str]] = None, ) -> TWEC: """Create a WEC object from linear hydrodynamic coefficients @@ -311,8 +311,8 @@ def from_bem( :py:func:`wecopttool.run_bem`, rather than running Capytaine directly, which outputs the results in the correct convention. The results can be saved - using :py:func:`wecopttool.write_netcdf`. - :py:func:`wecopttool.run_bem` also computes the inertia and + using :py:func:`wecopttool.write_netcdf`. + :py:func:`wecopttool.run_bem` also computes the inertia and hydrostatic stiffness which should be included in bem_data. Parameters @@ -323,7 +323,7 @@ def from_bem( corrected. Also includes inertia and hydrostatic stiffness. friction Linear friction, in addition to radiation damping, of size - :python:`(nodf, ndof)`. + :python:`(ndof, ndof)`. :python:`None` if included in :python:`bem_data` or to set to zero. f_add @@ -339,12 +339,12 @@ def from_bem( If :python:`None`: empty list :python:`[]`. min_damping Minimum damping level to ensure a stable system. - See :py:func:`wecopttool.check_linear_damping` for more details. + See :py:func:`wecopttool.check_radiation_damping` for more details. uniform_shift Boolean determining whether damping corrections shifts the damping values uniformly for all frequencies or only for frequencies below :python:`min_damping`. - See :py:func:`wecopttool.check_linear_damping` for more details. + See :py:func:`wecopttool.check_radiation_damping` for more details. dof_names Names of the different degrees of freedom (e.g. :python:`'Heave'`). @@ -354,7 +354,7 @@ def from_bem( See Also -------- run_bem, add_linear_friction, change_bem_convention, - write_netcdf, check_linear_damping + write_netcdf, check_radiation_damping """ if isinstance(bem_data, (str, Path)): bem_data = read_netcdf(bem_data) @@ -368,7 +368,7 @@ def from_bem( # check real part of damping diagonal > 0 if min_damping is not None: - hydro_data = check_linear_damping( + hydro_data = check_radiation_damping( hydro_data, min_damping, uniform_shift) # forces in the dynamics equations @@ -442,7 +442,7 @@ def from_floating_body( If :python:`None`: empty list :python:`[]`. min_damping Minimum damping level to ensure a stable system. - See :py:func:`wecopttool.check_linear_damping` for + See :py:func:`wecopttool.check_radiation_damping` for more details. wave_directions List of wave directions [degrees] to evaluate BEM at. @@ -483,6 +483,7 @@ def from_impedance( f_add: Optional[TIForceDict] = None, constraints: Optional[Iterable[Mapping]] = None, min_damping: Optional[float] = _default_min_damping, + uniform_shift: Optional[bool] = False, ) -> TWEC: """Create a WEC object from the intrinsic impedance and excitation coefficients. @@ -525,6 +526,11 @@ def from_impedance( Minimum damping level to ensure a stable system. See :py:func:`wecopttool.check_impedance` for more details. + uniform_shift + Boolean determining whether damping corrections shifts the damping + values uniformly for all frequencies or only for frequencies below + :python:`min_damping`. + See :py:func:`wecopttool.check_radiation_damping` for more details. Raises ------ @@ -688,18 +694,20 @@ def solve(self, may call >>> realization = 0 # realization index - >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt[realization]) - >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[realization]) + >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt,wave,nsubsteps) + >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt,wave,nsubsteps) See Also -------- wecopttool.waves, + wecopttool.core.wec.post_process, + wecopttool.core.pto.post_process, """ results = [] # x_wec scaling vector - if scale_x_wec == None: + if scale_x_wec is None: scale_x_wec = [1.0] * self.ndof elif isinstance(scale_x_wec, float) or isinstance(scale_x_wec, int): scale_x_wec = [scale_x_wec] * self.ndof @@ -714,13 +722,13 @@ def solve(self, # composite scaling vector scale = np.concatenate([scale_x_wec, scale_x_opt]) - + # decision variable initial guess if x_wec_0 is None: x_wec_0 = np.random.randn(self.nstate_wec) if x_opt_0 is None: x_opt_0 = np.random.randn(nstate_opt) - x0 = np.concatenate([x_wec_0, x_opt_0])*scale + x0 = np.concatenate([x_wec_0, x_opt_0])*scale # bounds if (bounds_wec is None) and (bounds_opt is None): @@ -748,11 +756,11 @@ def solve(self, for realization, wave in waves.groupby('realization'): _log.info("Solving pseudo-spectral control problem " - + f"for realization number {realization}.") - + + f"for realization number {realization}.") + # objective function sign = -1.0 if maximize else 1.0 - + def obj_fun_scaled(x): x_wec, x_opt = self.decompose_state(x/scale) return obj_fun(self, x_wec, x_opt, wave)*scale_obj*sign @@ -784,21 +792,6 @@ def scaled_resid_fun(x): if use_grad: eq_cons['jac'] = jacobian(scaled_resid_fun) constraints.append(eq_cons) - - # callback - if callback is None: - def callback_scipy(x): - x_wec, x_opt = self.decompose_state(x) - max_x_opt = np.nan if np.size(x_opt)==0 else np.max(np.abs(x_opt)) - _log.info("Scaled [max(x_wec), max(x_opt), obj_fun(x)]: " - + f"[{np.max(np.abs(x_wec)):.2e}, " - + f"{max_x_opt:.2e}, " - + f"{obj_fun_scaled(x):.2e}]") - else: - def callback_scipy(x): - x_s = x/scale - x_wec, x_opt = self.decompose_state(x_s) - return callback(self, x_wec, x_opt, wave) # optimization problem optim_options['disp'] = optim_options.get('disp', True) @@ -808,13 +801,12 @@ def callback_scipy(x): 'constraints': constraints, 'options': optim_options, 'bounds': bounds, - 'callback': callback_scipy, } if use_grad: problem['jac'] = grad(obj_fun_scaled) # minimize - optim_res = minimize(**problem) + optim_res = minimize_ipopt(**problem) msg = f'{optim_res.message} (Exit mode {optim_res.status})' if optim_res.status == 0: @@ -828,26 +820,28 @@ def callback_scipy(x): optim_res.x = optim_res.x / scale optim_res.fun = optim_res.fun / scale_obj optim_res.jac = optim_res.jac / scale_obj * scale - + results.append(optim_res) - + return results def post_process(self, - res: OptimizeResult, + wec: TWEC, + res: Union[OptimizeResult, Iterable], waves: Dataset, nsubsteps: Optional[int] = 1, - ) -> tuple[Dataset, Dataset]: - """Post-process the results from - :py:meth:`wecopttool.WEC.solve`. + ) -> tuple[list[Dataset], list[Dataset]]: + """Post-process the results from :py:meth:`wecopttool.WEC.solve`. Parameters ---------- + wec + :py:class:`wecopttool.WEC` object. + res + Results produced by :py:meth:`wecopttool.WEC.solve`. waves :py:class:`xarray.Dataset` with the structure and elements shown by :py:mod:`wecopttool.waves`. - res - Results produced by :py:func:`scipy.optimize.minimize`. nsubsteps Number of steps between the default (implied) time steps. A value of :python:`1` corresponds to the default step @@ -869,95 +863,107 @@ def post_process(self, obj_fun=pto.average_power, nstate_opt=2*nfreq+1) - To get the post-processed results for the - :py:class:`wecopttool.pto.PTO`, you may call + To get the post-processed results we may call - >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt) - >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt) - """ - create_time = f"{datetime.utcnow()}" - - omega_vals = np.concatenate([[0], waves.omega.values]) - freq_vals = np.concatenate([[0], waves.freq.values]) - period_vals = np.concatenate([[np.inf], 1/waves.freq.values]) - pos_attr = {'long_name': 'Position', 'units': 'm or rad'} - vel_attr = {'long_name': 'Velocity', 'units': 'm/s or rad/s'} - acc_attr = {'long_name': 'Acceleration', 'units': 'm/s^2 or rad/s^2'} - omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'} - freq_attr = {'long_name': 'Frequency', 'units': 'Hz'} - period_attr = {'long_name': 'Period', 'units': 's'} - time_attr = {'long_name': 'Time', 'units': 's'} - dof_attr = {'long_name': 'Degree of freedom'} - force_attr = {'long_name': 'Force or moment', 'units': 'N or Nm'} - wave_elev_attr = {'long_name': 'Wave elevation', 'units': 'm'} - x_wec, x_opt = self.decompose_state(res.x) - omega_coord = ("omega", omega_vals, omega_attr) - freq_coord = ("omega", freq_vals, freq_attr) - period_coord = ("omega", period_vals, period_attr) - dof_coord = ("influenced_dof", self.dof_names, dof_attr) - - # frequency domain - force_da_list = [] - for name, force in self.forces.items(): - force_td_tmp = force(self, x_wec, x_opt, waves) - force_fd = self.td_to_fd(force_td_tmp) - force_da = DataArray(data=force_fd, - dims=["omega", "influenced_dof"], - coords={ - 'omega': omega_coord, - 'freq': freq_coord, - 'period': period_coord, - 'influenced_dof': dof_coord}, - attrs=force_attr - ).expand_dims({'type': [name]}) - force_da_list.append(force_da) - - fd_forces = xr.concat(force_da_list, dim='type') - fd_forces.type.attrs['long_name'] = 'Type' - fd_forces.name = 'force' - fd_forces.attrs['long_name'] = 'Force' - - pos = self.vec_to_dofmat(x_wec) - pos_fd = real_to_complex(pos) - - vel = self.derivative_mat @ pos - vel_fd = real_to_complex(vel) - - acc = self.derivative2_mat @ pos - acc_fd = real_to_complex(acc) - - fd_state = Dataset( - data_vars={ - 'pos': (['omega', 'influenced_dof'], pos_fd, pos_attr), - 'vel': (['omega', 'influenced_dof'], vel_fd, vel_attr), - 'acc': (['omega', 'influenced_dof'], acc_fd, acc_attr)}, - coords={ - 'omega': omega_coord, - 'freq': freq_coord, - 'period': period_coord, - 'influenced_dof': dof_coord}, - attrs={"time_created_utc": create_time} - ) - - results_fd = xr.merge([fd_state, fd_forces, waves]) - results_fd = results_fd.transpose('omega', 'influenced_dof', 'type', - 'wave_direction') - results_fd = results_fd.fillna(0) - - # time domain - t_dat = self.time_nsubsteps(nsubsteps) - time = DataArray( - data=t_dat, name='time', dims='time', coords=[t_dat]) - results_td = results_fd.map(lambda x: time_results(x, time)) - - results_td['pos'].attrs = pos_attr - results_td['vel'].attrs = vel_attr - results_td['acc'].attrs = acc_attr - results_td['wave_elev'].attrs = wave_elev_attr - results_td['force'].attrs = force_attr - results_td['time'].attrs = time_attr - results_td.attrs['time_created_utc'] = create_time + >>> res_wec_fd, res_wec_td = wec.post_process(wec, res_opt,wave) + Note that :py:meth:`wecopttool.WEC.solve` method produces a list of + results objects (one for each phase realization). + """ + assert self == wec , ("The same wec object should be used to call " + + "post-process and be passed as an input.") + + def _postproc(res, waves, nsubsteps): + create_time = f"{datetime.utcnow()}" + + omega_vals = np.concatenate([[0], waves.omega.values]) + freq_vals = np.concatenate([[0], waves.freq.values]) + period_vals = np.concatenate([[np.inf], 1/waves.freq.values]) + pos_attr = {'long_name': 'Position', 'units': 'm or rad'} + vel_attr = {'long_name': 'Velocity', 'units': 'm/s or rad/s'} + acc_attr = {'long_name': 'Acceleration', 'units': 'm/s^2 or rad/s^2'} + omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'} + freq_attr = {'long_name': 'Frequency', 'units': 'Hz'} + period_attr = {'long_name': 'Period', 'units': 's'} + time_attr = {'long_name': 'Time', 'units': 's'} + dof_attr = {'long_name': 'Degree of freedom'} + force_attr = {'long_name': 'Force or moment', 'units': 'N or Nm'} + wave_elev_attr = {'long_name': 'Wave elevation', 'units': 'm'} + x_wec, x_opt = self.decompose_state(res.x) + omega_coord = ("omega", omega_vals, omega_attr) + freq_coord = ("omega", freq_vals, freq_attr) + period_coord = ("omega", period_vals, period_attr) + dof_coord = ("influenced_dof", self.dof_names, dof_attr) + + # frequency domain + force_da_list = [] + for name, force in self.forces.items(): + force_td_tmp = force(self, x_wec, x_opt, waves) + force_fd = self.td_to_fd(force_td_tmp) + force_da = DataArray(data=force_fd, + dims=["omega", "influenced_dof"], + coords={ + 'omega': omega_coord, + 'freq': freq_coord, + 'period': period_coord, + 'influenced_dof': dof_coord}, + attrs=force_attr + ).expand_dims({'type': [name]}) + force_da_list.append(force_da) + + fd_forces = xr.concat(force_da_list, dim='type') + fd_forces.type.attrs['long_name'] = 'Type' + fd_forces.name = 'force' + fd_forces.attrs['long_name'] = 'Force' + + pos = self.vec_to_dofmat(x_wec) + pos_fd = real_to_complex(pos) + + vel = self.derivative_mat @ pos + vel_fd = real_to_complex(vel) + + acc = self.derivative2_mat @ pos + acc_fd = real_to_complex(acc) + + fd_state = Dataset( + data_vars={ + 'pos': (['omega', 'influenced_dof'], pos_fd, pos_attr), + 'vel': (['omega', 'influenced_dof'], vel_fd, vel_attr), + 'acc': (['omega', 'influenced_dof'], acc_fd, acc_attr)}, + coords={ + 'omega': omega_coord, + 'freq': freq_coord, + 'period': period_coord, + 'influenced_dof': dof_coord}, + attrs={"time_created_utc": create_time} + ) + + results_fd = xr.merge([fd_state, fd_forces, waves]) + results_fd = results_fd.transpose('omega', 'influenced_dof', 'type', + 'wave_direction') + results_fd = results_fd.fillna(0) + + # time domain + t_dat = self.time_nsubsteps(nsubsteps) + time = DataArray( + data=t_dat, name='time', dims='time', coords=[t_dat]) + results_td = results_fd.map(lambda x: time_results(x, time)) + + results_td['pos'].attrs = pos_attr + results_td['vel'].attrs = vel_attr + results_td['acc'].attrs = acc_attr + results_td['wave_elev'].attrs = wave_elev_attr + results_td['force'].attrs = force_attr + results_td['time'].attrs = time_attr + results_td.attrs['time_created_utc'] = create_time + return results_fd, results_td + + results_fd = [] + results_td = [] + for idx, ires in enumerate(res): + ifd, itd = _postproc(ires, waves.sel(realization=idx), nsubsteps) + results_fd.append(ifd) + results_td.append(itd) return results_fd, results_td # properties @@ -1746,7 +1752,8 @@ def fd_to_td( td = tmat @ complex_to_real(fd, zero_freq) elif (f1 is None) and (nfreq is None): n = 2*(fd.shape[0]-1) - td = np.fft.irfft(fd/2, n=n, axis=0, norm='forward') + fd = np.concatenate((fd[:1, :], fd[1:-1, :]/2, fd[-1:, :])) + td = np.fft.irfft(fd, n=n, axis=0, norm='forward') else: raise ValueError( "Provide either both 'f1' and 'nfreq' or neither.") @@ -1782,10 +1789,10 @@ def td_to_fd( td= atleast_2d(td) n = td.shape[0] if fft: - fd = np.fft.rfft(td*2, n=n, axis=0, norm='forward') + fd = np.fft.rfft(td, n=n, axis=0, norm='forward') else: - fd = np.dot(dft(n, 'n')[:n//2+1, :], td*2) - fd = np.concatenate((fd[:1, :]/2, fd[1:-1, :], fd[-1:, :]/2)) + fd = np.dot(dft(n, 'n')[:n//2+1, :], td) + fd = np.concatenate((fd[:1, :], fd[1:-1, :]*2, fd[-1:, :])) if not zero_freq: fd = fd[1:, :] return fd @@ -1833,10 +1840,10 @@ def write_netcdf(fpath: Union[str, Path], data: Dataset) -> None: cpy.io.xarray.separate_complex_values(data).to_netcdf(fpath) -def check_linear_damping( +def check_radiation_damping( hydro_data: Dataset, min_damping: Optional[float] = 1e-6, - uniform_shift: Optional[bool] = True, + uniform_shift: Optional[bool] = False, ) -> Dataset: """Ensure that the linear hydrodynamics (friction + radiation damping) have positive damping. @@ -1858,7 +1865,7 @@ def check_linear_damping( damping for all frequencies. If :python:`False`, the damping correction is applied to :python:`radiation_damping` and only shifts the damping for frequencies with negative damping values. Default is - :python:`True`. + :python:`False`. """ hydro_data_new = hydro_data.copy(deep=True) radiation = hydro_data_new['radiation_damping'] @@ -1875,9 +1882,9 @@ def check_linear_damping( delta = min_damping-dmin _log.warning( f'Linear damping for DOF "{dof}" has negative or close ' + - 'to zero terms. Shifting up via linear friction of ' + + 'to zero terms. Shifting up radiation damping by ' + f'{delta.values} N/(m/s).') - hydro_data_new['friction'][idof, idof] = (ifriction + delta) + hydro_data_new['radiation_damping'][:, idof, idof] = (iradiation + delta) else: new_damping = iradiation.where( iradiation+ifriction>min_damping, other=min_damping) @@ -1885,7 +1892,8 @@ def check_linear_damping( if (new_damping==min_damping).any(): _log.warning( f'Linear damping for DOF "{dof}" has negative or close to ' + - 'zero terms. Shifting up damping terms to a minimum of ' + + 'zero terms. Shifting up damping terms ' + + f'{np.where(new_damping==min_damping)[0]} to a minimum of ' + f'{min_damping} N/(m/s)') hydro_data_new['radiation_damping'][:, idof, idof] = new_damping return hydro_data_new @@ -1894,6 +1902,7 @@ def check_linear_damping( def check_impedance( Zi: DataArray, min_damping: Optional[float] = 1e-6, + uniform_shift: Optional[bool] = False, ) -> DataArray: """Ensure that the real part of the impedance (resistive) is positive. @@ -1911,14 +1920,27 @@ def check_impedance( Zi_diag = np.diagonal(Zi,axis1=1,axis2=2) Zi_shifted = Zi.copy() for dof in range(Zi_diag.shape[1]): - dmin = np.min(np.real(Zi_diag[:, dof])) - if dmin < min_damping: - delta = min_damping - dmin - Zi_shifted[:, dof, dof] = Zi_diag[:, dof] \ - + np.abs(delta) - _log.warning( - f'Real part of impedance for {dof} has negative or close to ' + - f'zero terms. Shifting up by {delta:.2f}') + if uniform_shift: + dmin = np.min(np.real(Zi_diag[:, dof])) + if dmin < min_damping: + delta = min_damping - dmin + Zi_shifted[:, dof, dof] = Zi_diag[:, dof] \ + + np.abs(delta) + _log.warning( + f'Real part of impedance for {dof} has negative or close to ' + + f'zero terms. Shifting up by {delta:.2f}') + else: + points = np.where(np.real(Zi_diag[:, dof]) Dataset: :py:func:`wecopttool.add_linear_friction`. """ - Zi = (hydro_data['inertia_matrix'] \ + hydro_impedance = (hydro_data['inertia_matrix'] \ + hydro_data['added_mass'])*1j*hydro_data['omega'] \ + hydro_data['radiation_damping'] + hydro_data['friction'] \ + hydro_data['hydrostatic_stiffness']/1j/hydro_data['omega'] - return Zi.transpose('omega', 'radiating_dof', 'influenced_dof') + return hydro_impedance.transpose('omega', 'radiating_dof', 'influenced_dof') def atleast_2d(array: ArrayLike) -> ndarray: @@ -2521,7 +2543,7 @@ def set_fb_centers( - `rotation_center` is set to the center of mass """ valid_properties = ['center_of_mass', 'rotation_center'] - + for property in valid_properties: if not hasattr(fb, property): setattr(fb, property, None) @@ -2540,5 +2562,5 @@ def set_fb_centers( elif getattr(fb, property) is not None: _log.warning( f'{property} already defined as {getattr(fb, property)}.') - + return fb diff --git a/wecopttool/pto.py b/wecopttool/pto.py index e34b3b1e0..17841ec9c 100644 --- a/wecopttool/pto.py +++ b/wecopttool/pto.py @@ -656,10 +656,10 @@ def transduced_effort(self, def post_process(self, wec: TWEC, - res: OptimizeResult, + res: Union[OptimizeResult, list], waves: Optional[DataArray] = None, nsubsteps: Optional[int] = 1, - ) -> tuple[Dataset, Dataset]: + ) -> tuple[list[Dataset], list[Dataset]]: """Transform the results from optimization solution to a form that the user can work with directly. @@ -675,22 +675,21 @@ def post_process(self, To get the post-processed results for the :py:class:`wecopttool.pto.PTO`, you may call - >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt) - >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt) + >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[0],wave) For smoother plots, you can set :python:`nsubsteps` to a value greater than 1. >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt, nsubsteps=4) - >>> res_pto_td.power.plot() + >>> res_pto_td[0].power.plot() Parameters ---------- wec :py:class:`wecopttool.WEC` object. res - Results produced by :py:func:`scipy.optimize.minimize`. + Results produced by :py:meth:`wecopttool.WEC.solve`. waves :py:class:`xarray.Dataset` with the structure and elements shown by :py:mod:`wecopttool.waves`. @@ -702,115 +701,125 @@ def post_process(self, Returns ------- results_fd - :py:class:`xarray.Dataset` with frequency domain results. + list of :py:class:`xarray.Dataset` with frequency domain results. results_td - :py:class:`xarray.Dataset` with time domain results. + list of :py:class:`xarray.Dataset` with time domain results. """ - create_time = f"{datetime.utcnow()}" - - x_wec, x_opt = wec.decompose_state(res.x) - - # position - pos_td = self.position(wec, x_wec, x_opt, waves, nsubsteps) - pos_fd = wec.td_to_fd(pos_td[::nsubsteps]) - - # velocity - vel_td = self.velocity(wec, x_wec, x_opt, waves, nsubsteps) - vel_fd = wec.td_to_fd(vel_td[::nsubsteps]) - - # acceleration - acc_td = self.acceleration(wec, x_wec, x_opt, waves, nsubsteps) - acc_fd = wec.td_to_fd(acc_td[::nsubsteps]) - - # force - force_td = self.force(wec, x_wec, x_opt, waves, nsubsteps) - force_fd = wec.td_to_fd(force_td[::nsubsteps]) - - # power - elec_power_td = self.power(wec, x_wec, x_opt, waves, nsubsteps) - elec_power_fd = wec.td_to_fd(elec_power_td[::nsubsteps]) - - # mechanical power - mech_power_td = self.mechanical_power(wec, x_wec, x_opt, waves, - nsubsteps) - mech_power_fd = wec.td_to_fd(mech_power_td[::nsubsteps]) - - # stack mechanical and electrical power - power_names = ['mech','elec'] - power_fd = np.stack((mech_power_fd,elec_power_fd)) - power_td = np.stack((mech_power_td,elec_power_td)) - - pos_attr = {'long_name': 'Position', 'units': 'm or rad'} - vel_attr = {'long_name': 'Velocity', 'units': 'm/s or rad/s'} - acc_attr = {'long_name': 'Acceleration', - 'units': 'm/s^2 or rad/s^2'} - force_attr = {'long_name': 'Force or moment on WEC', - 'units': 'N or Nm'} - power_attr = {'long_name': 'Power', 'units': 'W'} - mech_power_attr = {'long_name': 'Mechanical power', 'units': 'W'} - omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'} - freq_attr = {'long_name': 'Frequency', 'units': 'Hz'} - period_attr = {'long_name': 'Period', 'units': 's'} - dof_attr = {'long_name': 'PTO degree of freedom'} - time_attr = {'long_name': 'Time', 'units': 's'} - type_attr = {'long_name': 'Power type'} - - t_dat = wec.time_nsubsteps(nsubsteps) - - results_fd = Dataset( - data_vars={ - 'pos': (['omega','dof'], pos_fd, pos_attr), - 'vel': (['omega','dof'], vel_fd, vel_attr), - 'acc': (['omega','dof'], acc_fd, acc_attr), - 'force': (['omega','dof'], force_fd, force_attr), - 'power': (['type','omega','dof'], power_fd, power_attr), - }, - coords={ - 'omega':('omega', wec.omega, omega_attr), - 'freq':('omega', wec.frequency, freq_attr), - 'period':('omega', wec.period, period_attr), - 'dof':('dof', self.names, dof_attr), - 'type':('type', power_names, power_attr)}, - attrs={"time_created_utc": create_time} - ) - - results_td = Dataset( - data_vars={ - 'pos': (['time','dof'], pos_td, pos_attr), - 'vel': (['time','dof'], vel_td, vel_attr), - 'acc': (['time','dof'], acc_td, acc_attr), - 'force': (['time','dof'], force_td, force_attr), - 'power': (['type','time','dof'], power_td, power_attr), - }, - coords={ - 'time':('time', t_dat, time_attr), - 'dof':('dof', self.names, dof_attr), - 'type':('type', power_names, power_attr)}, - attrs={"time_created_utc": create_time} - ) - - if self.impedance is not None: - #transduced flow and effort variables - q2_td, e2_td = self.power_variables(wec, x_wec, x_opt, - waves, nsubsteps) - q2_fd = wec.td_to_fd(q2_td[::nsubsteps]) - e2_fd = wec.td_to_fd(e2_td[::nsubsteps]) - - q2_attr = {'long_name': 'Transduced Flow', - 'units': 'A or m^3/s or rad/s or m/s'} - e2_attr = {'long_name': 'Transduced Effort', - 'units': 'V or N/m^2 or Nm or Ns'} - - results_td = results_td.assign({ - 'trans_flo': (['time','dof'], q2_td, q2_attr), - 'trans_eff': (['time','dof'], e2_td, e2_attr), - }) - results_fd = results_fd.assign({ - 'trans_flo': (['omega','dof'], q2_fd, q2_attr), - 'trans_eff': (['omega','dof'], e2_fd, e2_attr), - }) + def _postproc(wec, res, waves, nsubsteps): + + create_time = f"{datetime.utcnow()}" + + x_wec, x_opt = wec.decompose_state(res.x) + + # position + pos_td = self.position(wec, x_wec, x_opt, waves, nsubsteps) + pos_fd = wec.td_to_fd(pos_td[::nsubsteps]) + + # velocity + vel_td = self.velocity(wec, x_wec, x_opt, waves, nsubsteps) + vel_fd = wec.td_to_fd(vel_td[::nsubsteps]) + + # acceleration + acc_td = self.acceleration(wec, x_wec, x_opt, waves, nsubsteps) + acc_fd = wec.td_to_fd(acc_td[::nsubsteps]) + + # force + force_td = self.force(wec, x_wec, x_opt, waves, nsubsteps) + force_fd = wec.td_to_fd(force_td[::nsubsteps]) + + # power + elec_power_td = self.power(wec, x_wec, x_opt, waves, nsubsteps) + elec_power_fd = wec.td_to_fd(elec_power_td[::nsubsteps]) + + # mechanical power + mech_power_td = self.mechanical_power(wec, x_wec, x_opt, waves, + nsubsteps) + mech_power_fd = wec.td_to_fd(mech_power_td[::nsubsteps]) + + # stack mechanical and electrical power + power_names = ['mech','elec'] + power_fd = np.stack((mech_power_fd,elec_power_fd)) + power_td = np.stack((mech_power_td,elec_power_td)) + + pos_attr = {'long_name': 'Position', 'units': 'm or rad'} + vel_attr = {'long_name': 'Velocity', 'units': 'm/s or rad/s'} + acc_attr = {'long_name': 'Acceleration', + 'units': 'm/s^2 or rad/s^2'} + force_attr = {'long_name': 'Force or moment on WEC', + 'units': 'N or Nm'} + power_attr = {'long_name': 'Power', 'units': 'W'} + mech_power_attr = {'long_name': 'Mechanical power', 'units': 'W'} + omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'} + freq_attr = {'long_name': 'Frequency', 'units': 'Hz'} + period_attr = {'long_name': 'Period', 'units': 's'} + dof_attr = {'long_name': 'PTO degree of freedom'} + time_attr = {'long_name': 'Time', 'units': 's'} + type_attr = {'long_name': 'Power type'} + + t_dat = wec.time_nsubsteps(nsubsteps) + + results_fd = Dataset( + data_vars={ + 'pos': (['omega','dof'], pos_fd, pos_attr), + 'vel': (['omega','dof'], vel_fd, vel_attr), + 'acc': (['omega','dof'], acc_fd, acc_attr), + 'force': (['omega','dof'], force_fd, force_attr), + 'power': (['type','omega','dof'], power_fd, power_attr), + }, + coords={ + 'omega':('omega', wec.omega, omega_attr), + 'freq':('omega', wec.frequency, freq_attr), + 'period':('omega', wec.period, period_attr), + 'dof':('dof', self.names, dof_attr), + 'type':('type', power_names, power_attr)}, + attrs={"time_created_utc": create_time} + ) + results_td = Dataset( + data_vars={ + 'pos': (['time','dof'], pos_td, pos_attr), + 'vel': (['time','dof'], vel_td, vel_attr), + 'acc': (['time','dof'], acc_td, acc_attr), + 'force': (['time','dof'], force_td, force_attr), + 'power': (['type','time','dof'], power_td, power_attr), + }, + coords={ + 'time':('time', t_dat, time_attr), + 'dof':('dof', self.names, dof_attr), + 'type':('type', power_names, power_attr)}, + attrs={"time_created_utc": create_time} + ) + if self.impedance is not None: + #transduced flow and effort variables + q2_td, e2_td = self.power_variables(wec, x_wec, x_opt, + waves, nsubsteps) + q2_fd = wec.td_to_fd(q2_td[::nsubsteps]) + e2_fd = wec.td_to_fd(e2_td[::nsubsteps]) + + q2_attr = {'long_name': 'Transduced Flow', + 'units': 'A or m^3/s or rad/s or m/s'} + e2_attr = {'long_name': 'Transduced Effort', + 'units': 'V or N/m^2 or Nm or Ns'} + + results_td = results_td.assign({ + 'trans_flo': (['time','dof'], q2_td, q2_attr), + 'trans_eff': (['time','dof'], e2_td, e2_attr), + }) + results_fd = results_fd.assign({ + 'trans_flo': (['omega','dof'], q2_fd, q2_attr), + 'trans_eff': (['omega','dof'], e2_fd, e2_attr), + }) + + + return results_fd, results_td + + results_fd = [] + results_td = [] + for idx, ires in enumerate(res): + ifd, itd = _postproc(wec, ires, waves.sel(realization=idx), nsubsteps) + results_fd.append(ifd) + results_td.append(itd) return results_fd, results_td diff --git a/wecopttool/utilities.py b/wecopttool/utilities.py new file mode 100644 index 000000000..d76731704 --- /dev/null +++ b/wecopttool/utilities.py @@ -0,0 +1,369 @@ +"""Functions that are useful for WEC analysis and design. +""" + + +from __future__ import annotations + + +__all__ = [ + "plot_hydrodynamic_coefficients", + "plot_bode_impedance", + "calculate_power_flows", + "plot_power_flow", +] + + +from typing import Optional, Union +import logging +from pathlib import Path + +import autograd.numpy as np +from autograd.numpy import ndarray + +from xarray import DataArray +from numpy.typing import ArrayLike +# from autograd.numpy import ndarray +from xarray import DataArray, concat +import matplotlib.pyplot as plt +from matplotlib.figure import Figure +from matplotlib.axes import Axes + +from matplotlib.sankey import Sankey + + + +# logger +_log = logging.getLogger(__name__) + + + + +def plot_hydrodynamic_coefficients(bem_data, + wave_dir: Optional[float] = 0.0 + )-> list(tuple(Figure, Axes)): + """Plots hydrodynamic coefficients (added mass, radiation damping, + and wave excitation) based on BEM data. + + + Parameters + ---------- + bem_data + Linear hydrodynamic coefficients obtained using the boundary + element method (BEM) code Capytaine, with sign convention + corrected. + wave_dir + Wave direction(s) to plot the + + """ + + bem_data = bem_data.sel(wave_direction = wave_dir, method='nearest') + radiating_dofs = bem_data.radiating_dof.values + influenced_dofs = bem_data.influenced_dof.values + + # plots + fig_am, ax_am = plt.subplots( + len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs),3*len(influenced_dofs)), + squeeze=False + ) + fig_rd, ax_rd = plt.subplots( + len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs), 3*len(influenced_dofs)), + squeeze=False + ) + fig_ex, ax_ex = plt.subplots( + len(influenced_dofs), + 1, + tight_layout=True, + sharex=True, + figsize=(3, 3*len(radiating_dofs)), + squeeze=False + ) + [ax.grid(True) for axs in (ax_am, ax_rd, ax_ex) for ax in axs.flatten()] + # plot titles + fig_am.suptitle('Added Mass Coefficients', fontweight='bold') + fig_rd.suptitle('Radiation Damping Coefficients', fontweight='bold') + fig_ex.suptitle('Wave Excitation Coefficients', fontweight='bold') + + sp_idx = 0 + for i, rdof in enumerate(radiating_dofs): + for j, idof in enumerate(influenced_dofs): + sp_idx += 1 + if i == 0: + np.abs(bem_data.diffraction_force.sel(influenced_dof=idof)).plot( + ax=ax_ex[j,0], linestyle='dashed', label='Diffraction') + np.abs(bem_data.Froude_Krylov_force.sel(influenced_dof=idof)).plot( + ax=ax_ex[j,0], linestyle='dashdot', label='Froude-Krylov') + ex_handles, ex_labels = ax_ex[j,0].get_legend_handles_labels() + ax_ex[j,0].set_title(f'{idof}') + ax_ex[j,0].set_xlabel('') + ax_ex[j,0].set_ylabel('') + if j <= i: + bem_data.added_mass.sel( + radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_am[i, j]) + bem_data.radiation_damping.sel( + radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_rd[i, j]) + if i == len(radiating_dofs)-1: + ax_am[i, j].set_xlabel(f'$\omega$', fontsize=10) + ax_rd[i, j].set_xlabel(f'$\omega$', fontsize=10) + ax_ex[j, 0].set_xlabel(f'$\omega$', fontsize=10) + else: + ax_am[i, j].set_xlabel('') + ax_rd[i, j].set_xlabel('') + if j == 0: + ax_am[i, j].set_ylabel(f'{rdof}', fontsize=10) + ax_rd[i, j].set_ylabel(f'{rdof}', fontsize=10) + else: + ax_am[i, j].set_ylabel('') + ax_rd[i, j].set_ylabel('') + if j == i: + ax_am[i, j].set_title(f'{idof}', fontsize=10) + ax_rd[i, j].set_title(f'{idof}', fontsize=10) + else: + ax_am[i, j].set_title('') + ax_rd[i, j].set_title('') + else: + fig_am.delaxes(ax_am[i, j]) + fig_rd.delaxes(ax_rd[i, j]) + fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False) + return [(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)] + +def plot_bode_impedance(impedance: DataArray, + title: Optional[str]= '', + #plot_natural_freq: Optional[bool] = False, +)-> tuple(Figure, Axes): + """Plot Bode graph from wecoptool impedance data array. + + Parameters + ---------- + impedance: DataArray + Complex impedance matrix produced by for example by + :py:func:`wecopttool.hydrodynamic_impedance`. + Dimensions: omega, radiating_dofs, influenced_dofs + title: String + Title string to be displayed in the plot. + + """ + radiating_dofs = impedance.radiating_dof.values + influenced_dofs = impedance.influenced_dof.values + mag = 20.0 * np.log10(np.abs(impedance)) + phase = np.rad2deg(np.unwrap(np.angle(impedance))) + freq = impedance.omega.values/2/np.pi + fig, axes = plt.subplots( + 2*len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs), 3*len(influenced_dofs)), + squeeze=False + ) + fig.suptitle(title + ' Bode Plots', fontweight='bold') + + sp_idx = 0 + for i, rdof in enumerate(radiating_dofs): + for j, idof in enumerate(influenced_dofs): + sp_idx += 1 + axes[2*i, j].semilogx(freq, mag[:, i, j]) # Bode magnitude plot + axes[2*i+1, j].semilogx(freq, phase[:, i, j]) # Bode phase plot + axes[2*i, j].grid(True, which = 'both') + axes[2*i+1, j].grid(True, which = 'both') + if i == len(radiating_dofs)-1: + axes[2*i+1, j].set_xlabel(f'Frequency (Hz)', fontsize=10) + else: + axes[i, j].set_xlabel('') + if j == 0: + axes[2*i, j].set_ylabel(f'{rdof} \n Mag. (dB)', fontsize=10) + axes[2*i+1, j].set_ylabel(f'Phase. (deg)', fontsize=10) + else: + axes[i, j].set_ylabel('') + if i == 0: + axes[i, j].set_title(f'{idof}', fontsize=10) + else: + axes[i, j].set_title('') + return fig, axes + + + + +def calculate_power_flows(wec, + pto, + results, + waves, + intrinsic_impedance)-> dict[str, float]: + """Calculate power flows into a :py:class:`wecopttool.WEC` + and through a :py:class:`wecopttool.pto.PTO` based on the results + of :py:meth:`wecopttool.WEC.solve` for a single wave realization. + + Parameters + ---------- + wec + WEC object of :py:class:`wecopttool.WEC` + pto + PTO object of :py:class:`wecopttool.pto.PTO` + results + Results produced by :py:func:`scipy.optimize.minimize` for a single wave + realization. + waves + :py:class:`xarray.Dataset` with the structure and elements + shown by :py:mod:`wecopttool.waves`. + intrinsic_impedance: DataArray + Complex intrinsic impedance matrix produced by + :py:func:`wecopttool.hydrodynamic_impedance`. + Dimensions: omega, radiating_dofs, influenced_dofs + + + """ + wec_fdom, _ = wec.post_process(wec, results, waves) + x_wec, x_opt = wec.decompose_state(results[0].x) + + #power quntities from solver + P_mech = pto.mechanical_average_power(wec, x_wec, x_opt, waves) + P_elec = pto.average_power(wec, x_wec, x_opt, waves) + + #compute analytical power flows + Fex_FD = wec_fdom[0].force.sel(type=['Froude_Krylov', 'diffraction']).sum('type') + Rad_res = np.real(intrinsic_impedance.squeeze()) + Vel_FD = wec_fdom[0].vel + + P_max, P_e, P_r = [], [], [] + + #This solution requires radiation resistance matrix Rad_res to be invertible + # TODO In the future we might want to add an entirely unconstrained solve + # for optimized mechanical power + + for om in Rad_res.omega.values: + #use frequency vector from intrinsic impedance (no zero freq) + #Eq. 6.69 + #Dofs are row vector, which is transposed in standard convention + Fe_FD_t = np.atleast_2d(Fex_FD.sel(omega = om)) + Fe_FD = np.transpose(Fe_FD_t) + R_inv = np.linalg.inv(np.atleast_2d(Rad_res.sel(omega= om))) + P_max.append((1/8)*(Fe_FD_t@R_inv)@np.conj(Fe_FD)) + #Eq.6.57 + U_FD_t = np.atleast_2d(Vel_FD.sel(omega = om)) + U_FD = np.transpose(U_FD_t) + R = np.atleast_2d(Rad_res.sel(omega= om)) + P_r.append((1/2)*(U_FD_t@R)@np.conj(U_FD)) + #Eq. 6.56 (replaced pinv(Fe)*U with U'*conj(Fe) + # as suggested in subsequent paragraph) + P_e.append((1/4)*(Fe_FD_t@np.conj(U_FD) + U_FD_t@np.conj(Fe_FD))) + + power_flows = { + 'Optimal Excitation' : -2* np.sum(np.real(P_max)),#eq 6.68 + 'Radiated': -1*np.sum(np.real(P_r)), + 'Actual Excitation': -1*np.sum(np.real(P_e)), + 'Electrical (solver)': P_elec, + 'Mechanical (solver)': P_mech, + } + + power_flows['Absorbed'] = ( + power_flows['Actual Excitation'] + - power_flows['Radiated'] + ) + power_flows['Unused Potential'] = ( + power_flows['Optimal Excitation'] + - power_flows['Actual Excitation'] + ) + power_flows['PTO Loss'] = ( + power_flows['Mechanical (solver)'] + - power_flows['Electrical (solver)'] + ) + return power_flows + + +def plot_power_flow(power_flows: dict[str, float])-> tuple(Figure, Axes): + """Plot power flow through a WEC as Sankey diagram. + + Parameters + ---------- + power_flows: dictionary + Power flow dictionary produced by for example by + :py:func:`wecopttool.utilities.calculate_power_flows`. + Required keys: 'Optimal Excitation', 'Radiated', 'Actual Excitation' + 'Electrical (solver)', 'Mechanical (solver)', + 'Absorbed', 'Unused Potential', 'PTO Loss' + + + """ + # fig = plt.figure(figsize = [8,4]) + # ax = fig.add_subplot(1, 1, 1,) + fig, ax = plt.subplots(1, 1, + tight_layout=True, + figsize=(8, 4), + ) + + # plt.viridis() + sankey = Sankey(ax=ax, + scale= -1/power_flows['Optimal Excitation'], + offset= 0, + format = '%.1f', + shoulder = 0.02, + tolerance=-1e-03*power_flows['Optimal Excitation'], + unit = 'W' + ) + + sankey.add(flows=[-1*power_flows['Optimal Excitation'], + power_flows['Unused Potential'], + power_flows['Actual Excitation']], + labels = ['Optimal Excitation', + 'Unused Potential ', + 'Excited'], + orientations=[0, -1, -0],#arrow directions, + pathlengths = [0.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.253935, 0.265254, 0.529983, 1.0) #viridis(0.2) + ) + + sankey.add(flows=[ + -1*(power_flows['Absorbed'] + power_flows['Radiated']), + power_flows['Radiated'], + power_flows['Absorbed'], + ], + labels = ['Excited', + 'Radiated', + ''], + prior= (0), + connect=(2,0), + orientations=[0, -1, -0],#arrow directions, + pathlengths = [0.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.127568, 0.566949, 0.550556, 1.0) #viridis(0.5) + ) + + sankey.add(flows=[-1*(power_flows['Mechanical (solver)']), + power_flows['PTO Loss'], + power_flows['Electrical (solver)'], + ], + labels = ['Mechanical', + 'PTO-Loss' , + 'Electrical'], + prior= (1), + connect=(2,0), + orientations=[0, -1, -0],#arrow directions, + pathlengths = [.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.741388, 0.873449, 0.149561, 1.0) #viridis(0.9) + ) + + + diagrams = sankey.finish() + for diagram in diagrams: + for text in diagram.texts: + text.set_fontsize(10) + #remove text label from last entries + for diagram in diagrams[0:2]: + diagram.texts[2].set_text('') + + plt.axis("off") + # plt.show() + + return fig, ax diff --git a/wecopttool/waves.py b/wecopttool/waves.py index 97ddc1a29..0f33f88ef 100644 --- a/wecopttool/waves.py +++ b/wecopttool/waves.py @@ -5,13 +5,14 @@ :python:`wecopttool`. It also provides functions for creating common types of waves such as regular waves and irregular waves. -The data structure is a 2D complex :py:module:xarray.DataArray +The data structure is a 2D complex :py:class:`xarray.DataArray` containing the complex amplitude. The 2D coordinates are: wave angular frequency :python:`omega` (rad/s) and direction :python:`wave_direction` (rad). This module uses wave spectrum data in the -:py:class:`wavespectra.SpecArray` format. +:py:class:`wavespectra.SpecArray` format, but does not require that you +use :py:class:`wavespectra.SpecArray` objects. """ @@ -203,8 +204,8 @@ def long_crested_wave( """Create a complex frequency-domain wave elevation from an omnidirectional spectrum. - The omnidirectional spectrum is in the - :py:class:`wavespectra.SpecArray` format. + The spectrum is a :py:class:`xarray.DataArray` in the format used + by :py:class:`wavespectra.SpecArray`. .. note:: The frequencies must be evenly-spaced with spacing equal to the first frequency. This is not always the case when @@ -245,8 +246,8 @@ def irregular_wave(efth: DataArray, seed: Optional[float] = None,) -> DataArray: """Create a complex frequency-domain wave elevation from a spectrum. - The omnidirectional spectrum is in the - :py:class:`wavespectra.SpecArray` format. + The spectrum is a :py:class:`xarray.DataArray` in the format used + by :py:class:`wavespectra.SpecArray`. .. note:: The frequencies must be evenly-spaced with spacing equal to the first frequency. This is not always the case when @@ -308,8 +309,8 @@ def omnidirectional_spectrum( spectrum_func: Callable, spectrum_name: str = '', ) -> DataArray: - """Create the dataset for an omnidirectional wave spectrum in the - :py:class:`wavespectra.SpecArray` format. + """Create the :py:class:`xarray.DataArray` for an omnidirectional + wave spectrum in the :py:class:`wavespectra.SpecArray` format. Examples -------- @@ -378,8 +379,8 @@ def spectrum( spectrum_name: str = '', spread_name: str = '', ) -> DataArray: - """Create the dataset for an irregular wave in the - :py:class:`wavespectra.SpecArray` format. + """Create the :py:class:`xarray.DataArray` for an irregular wave + in the :py:class:`wavespectra.SpecArray` format. Examples --------