sandialabs / wecopttool

WEC Design Optimization Toolbox

WecOptTool

The Wave Energy Converter Design Optimization Toolbox (WecOptTool) allows users to perform wave energy converter (WEC) device design optimization studies with constrained optimal control.

NOTE: If you are looking for the WecOptTool code used in previous published work (MATLAB version) please see WecOptTool-MATLAB.

Project Information

Refer to WecOptTool documentation for more information, including project overview, tutorials, theory, and API documentation.

Getting started

If you are brand new to Python and/or want detailed installation instructions, click here.

WecOptTool requires Python >= 3.8. Python 3.10 & 3.11 are supported. It is strongly recommended you create a dedicated virtual environment (e.g., using conda, mamba, venv, etc.) before installing WecOptTool.

From your dedicated environment, you can install WecOptTool via conda, pip, or mamba:

Option 1 - using Conda:

conda install -c conda-forge wecopttool

Option 2 - using pip (requires Fortran compilers on your system):

pip install wecopttool

Option 3 - using Mamba:

mamba install wecopttool

Geometry module and tutorials

To use our geometry examples, including for running the tutorials, you will need to install some additional dependencies. For the tutorials you will also need to install jupyter.

pip install wecopttool[geometry] jupyter

or on a Mac (Zsh shell)

pip install wecopttool\[geometry] jupyter

Tutorials

The tutorials can be found in the examples directory and are written as Jupyter Notebooks. To run the tutorials, first download the notebook files and then, from the directory containing the notebooks, run jupyter notebook. Using git to obtain the notebooks this can be done by running

git clone https://github.com/sandialabs/WecOptTool.git
cd WecOptTool/examples
jupyter notebook

Getting help

To report bugs, use WecOptTool's issues page. For general discussion, use WecOptTool's discussion page

Contributing

If you are interested in contributing to WecOptTool, see our contribution guidelines.

wecopttool's Issues

bounds_wec and bounds_opt separately

With #71, the user can set bounds on the decision variable. We may want to only allow the user to set bounds on x_opt, not x_wec. Not entirely sure about this... let's hash it out here.

Intuitive hydrostatics

With #49, we have a conistent way of handling's hydrostatics, buoyancy. However, want to consider:

• is this as intuitive as possible
• can we use capytaine's upcoming hydrostatics module (capytaine/capytaine#106)?

Initialize WEC with impedance rather than a Capytaine floating body.

If the user already has a linear model (impedance) they shouldn't have to provide a mesh and DOFs (capytaine.floatingbody) or solve the BEM. This would also allow for easier implementation of non-linear hydrodynamics: could even set the linear impedance, or parts of it to zero, and specify the non-linear functions through f_add. This would be a major refactoring/API change.

BUG: install does not give "jupiter-notebook"

Describe the bug
There are two issues with our installation configuration:

1. The default installation does not include jupyter-notebook
2. The options ([dev] and [doc]) contain notebook, but should include jupyter (https://docs.jupyter.org/en/latest/install/notebook-classic.html#alternative-for-experienced-python-users-installing-jupyter-with-pip)

https://github.com/SNL-WaterPower/WecOptTool/blob/6bd7c8f8561a9bfdaf2b0b39f3b6921a1b1d8562/setup.cfg#L20-L50

To Reproduce
Steps to reproduce the behavior: follow install instructions on dREADME.md`:

https://github.com/SNL-WaterPower/WecOptTool/blob/6bd7c8f8561a9bfdaf2b0b39f3b6921a1b1d8562/README.md?plain=1#L24-L28

Expected behavior
Installation includes jupyter-notebook

Observed behavior

$ conda list note
# packages in environment at /Users/rcoe/anaconda3/envs/demo:
#
# Name                    Version                   Build  Channel
notebook                  6.4.8                    pypi_0    pypi

$ conda list jupyter
# packages in environment at /Users/rcoe/anaconda3/envs/demo:
#
# Name                    Version                   Build  Channel
jupyter-client            7.1.2                    pypi_0    pypi
jupyter-core              4.9.1                    pypi_0    pypi
jupyterlab-pygments       0.1.2                    pypi_0    pypi

System:

• OS: macOs
• Python version: 3.9.10
• WecOptTool version: 1.0.1

Feature request: use previous solution for initial guess

Feature description.
When running design optimization studies, it is often the case that you will be calling WecOptTool.solve for a number of different design configurations (e.g., slightly different geometries). If the optimization algorithm you're using executes these calls serially, you can use the output of one wec.solve to set the initial guess (x_wec_0 and x_opt_0) for the decision vector used in the subsequent wec.solve call. While the solution of the pseudo-spectral problem will change based on the wec, if the wec's are similar, as is often the case, the initial guess will be close to the optimal solution of x and thus reduce computational time.

Describe the solution you'd like
Probably no need to alter the source for this; just show how this can be done via examples

Feature request: Allow wave directions to be a subset of the BEM directions.

Feature description.
Allows the user to define multi-directional waves as long as the incident wave directions are a subset of the wave directions in the BEM solution

Issue addressed.
Would also solve #80.

Validate advanced features

The RM3 example currently runs but needs to be validated. It uses many of the advanced features not included in the WaveBot tutorials or validation case. These include:

• irregular waves
• multi-directional waves
• multi-body WEC
• multiple WEC degrees of freedom
• PTO degree of freedom is function of more than one WEC DOF

Feature request: Construct PTO Jacobian for rigid body motions

Feature description.
If the WEC DOFs are translational and rotational rigid body motions we can help the user construct the Jacobian matrix by asking them the type of PTO (translation, rotation) & the rotation axis for rotational. @dtgaebe has derived these equations. This would fit nicely in the PTO module.

Issue addressed.
Make it easier for users to use our PTO classes. Currently they have to calculate the Jacobian matrix and provide that as input, but we can easily construct this matrix for them.

Describe the solution you'd like
A function in the PTO module that users can optionally use.

Describe alternatives you've considered
Currently creating the Jacobian falls completely on the user.

Interest in leading this feature development?
@dtgaebe is interested in leading this.

Additional information

PTO <=> WEC degrees of freedom

I think the conversion between PTO and WEC DoF is wrong. Currently we're using a kinematics matrix to go in one direction, and the transpose of that matrix to go in the other. If it was merely a change of coordinates, the two matrices would be the inverse of each other, not transpose. However it is not simply a change of coordinate, the number of WEC DoF doesn't need to be the same as the PTO DoF (and even if a square matrix doesn't have to be nonsingular). Seems to me that two different transformations are needed and one cannot be inferred from the other.

1. From WEC DoF (solved in the dynamics) to PTO DoF (needed to calculate PTO forces)
2. From PTO forces to forces on the WEC DoF (needed to solve the dynamics)

BUG: multiple constraints, only one enforced...?

alternative optimization package?

We are currently using scipy.optimize.minimize.

We have previously played with using IPOPT for the pseudo-spectral solution. From what I can tell, this is developed using C++ by COIN-OR; there are a couple of Python wrappers.

Verify wave directions are close to direction in BEM

Feature description.
When doing wec.read_bem() it does not check for wave direction. It should check for wave direction in a similar way like checking for mass

Refactor core.py, properters, etc.

Our current system for setting properties and updating the wec class does what we want, but it is not easy to read/edit.

Probably tackle #84 with this.

Move to capytaine 1.2.1

I help make some minor updates to capytaine that are part of v1.2.1. In particular, there was a bug about the naming of the mass matrix. We should use this latest version of capytaine in WecOptTool.

changelog: https://github.com/mancellin/capytaine/blob/master/docs/changelog.rst

@cmichelenstrofer - as a general mode of operating, it may often sense to make alterations directly to capytaine instead of in WecOptTool - I definitely support this

expose callbacks argument to scipy.optimize.minimize

related to #75, we want to expose the callbacks argument for scipy.optimize.minimize

Natural frequency in multiple DOFs

Similarly to #69, we need to verify the functionality of our natural_frequency method and function in multiple degrees of freedom.

https://github.com/SNL-WaterPower/WecOptTool/blob/7d9bbc70f7cac23b468d5b073d6da9727c68be9e/wecopttool/core.py#L1105-L1127

Parallel Computation

capytaine/capytaine#65 (comment)

@ryancoe I was able to see the use of all of my cores when running the capytaine example.

Using the current WecOptTool main repository example_WaveBot_oneDof.py I am not sure the problem is difficult enough to view the issue. Can you tell what specifically you were running/ testing in WecOptTool to test this?

Test fixtures behavior

In the test_linear_pi_pto test: overwriting the values of the wec fixture had no effect, the original f_add is still used (PS PTO force, instead of Linear PI PTO force). I had to create a new _wec fixture that never had f_add to begin with. This might be an issue with some of the other tests (as code will still run, just with the original PTO force). Need to understand what is going on in this case and how fixtures are supposed to behave.

ensure real part of diagonals of impedance are positive

Due to irregular frequencies and perhaps some other numerical issues, it is possible to get negative values for the radiation damping. This can be valid on the off-diagonals (i.e., the terms that describe interaction between modes of motion and/or different bodies), but it cannot be valid for the diagonals. When negative damping occurs on the diagonals, you no longer have a passive system.

We currently have a default value for friction damping that is added to the diagonals:

https://github.com/SNL-WaterPower/WecOptTool/blob/3933c5ff9dd7320e2dea14a29233d40d32d4d880/WecOptTool/core.py#L75-L79

We should add an additional step that checks whether the diagonals of the impedance matrix have positive real parts, and sets those which are less than zero to zero. It may be good to return a warning that this is happening also.

BUG: Funky read_bem error messages

Describe the bug
Changes to mesh (and other things) can cause the usage of read_bem to throw an unhelpful error message.

To Reproduce
Steps to reproduce the behavior:

1. Run tutorial_3_optimization.ipynb
2. Adjust the mesh_size_factor
3. Re-run the notebook
4. This will throw the assertion error from the following line

https://github.com/SNL-WaterPower/WecOptTool/blob/19a523f65165a6ef7544613ab1c8336db0dd96d1/wecopttool/core.py#L391

Expected behavior
Better error message (e.g., "Did you change your mesh?")

System:

• OS: macOS
• Python version: 3.9.7
• WecOptTool version: 1.0.0
• Any other relevant package version: N/A

Additional information
Suggested solution: wrap the assertions in read_bem within a try-catch environment and add a more helpful error message when any of them fail

Validate irregular waves (directional and multidirectional)

Create validation cases that use irregular waves, both unidirectional and multidirectional. Compare to linear solution.

Power Limit in multiple DOFs

We need to verify that the power limit function works for multiple DOFs.

Explicitly consider buoyancy vs. gravity

Currently, we assume equilibrium for the gravity and buoyancy of the body. This can be a bit unintuitive. We want to alter the residual equation to include gravity and buoyancy. In the case that these are not balanced (e.g., a tether or mooring system provides some additional pretension), the user will need to include a f_add that balances the gravity and buoyancy forces.

https://github.com/SNL-WaterPower/WecOptTool/blob/51c3aea47098ab1a7cf6c2a8f9801ff60de5645e/wecopttool/core.py#L786

f_b = f_b0 + kx, where k = 0 for a submerged body

We considered altering the impedance, but gravity is not position (state) dependent.

BUG: P and PI PTOs for multiple DOFs

This might be broken for multiple DOFs due to wrong reshaping.

[File ~/Documents/WecOptTool/WecOptTool/wecopttool/pto.py:333, in ProportionalPTO.force(self, wec, x_wec, x_opt, nsubsteps)
    ]()[330]()[ def force(self, wec: WEC, x_wec: npt.ArrayLike, x_opt: npt.ArrayLike,
    ]()[331]()[           nsubsteps: int = 1) -> np.ndarray:
    ]()[332]()[     vel_td = self.velocity(wec, x_wec, x_opt, nsubsteps)
--> ]()[333]()[     force_td = np.reshape(x_opt, [-1, 1]) * vel_td
    ]()[334]()     return force_td

ValueError: operands could not be broadcast together with shapes (2,1) (101,2)

Unit Tests

The tests currently just test that that the code runs. Now that we have an initial release we should write tests that check that new commits do not inadvertently change the code outputs.

• use assert to verify new versions give same results as previous ones.

Refactor PTO module

• Fewer classes
• Better names
• More intuitive interface

Validate multiple WEC and PTO degrees of freedom

After #52 is resolved validate that the code produces correct results when the WEC has multiple DoF and the (one or more) PTO depend on a combination of these. This could be the WaveBot with 3 DoF and 2 PTOs. The results should be compared to linear theory to validate the solution to #52.

• validate code with linear theory for multiple DoF
• include clean example in the examples folder

Jinja ->v3.1

Jinja released v3.1.0 earlier today and it seems to create a build error for us related to including the Jupiter notebooks in the docs. The error message is module 'jinja2.utils' has no attribute 'escape'

2022-03-24T18:58:04.7554931Z ##[group]Run cd docs
2022-03-24T18:58:04.7555234Z �[36;1mcd docs�[0m
2022-03-24T18:58:04.7555488Z �[36;1mmake clean html�[0m
2022-03-24T18:58:04.7555746Z �[36;1mmake clean linkcheck�[0m
2022-03-24T18:58:04.7556005Z �[36;1mcd ../�[0m
2022-03-24T18:58:04.7610494Z shell: /usr/bin/bash -l {0}
2022-03-24T18:58:04.7610765Z env:
2022-03-24T18:58:04.7611049Z   CONDA_PKGS_DIR: /home/runner/conda_pkgs_dir
2022-03-24T18:58:04.7611326Z ##[endgroup]
2022-03-24T18:58:05.1157908Z python "source"/make_theory_animations.py
2022-03-24T18:58:17.4931774Z Running Sphinx v4.4.0
2022-03-24T18:58:20.3663321Z Copy example notebooks into docs/_examples
2022-03-24T18:58:20.3664236Z checking bibtex cache... out of date
2022-03-24T18:58:20.3694141Z parsing bibtex file /home/runner/work/WecOptTool/WecOptTool/docs/source/wecopttool_refs.bib... parsed 5 entries
2022-03-24T18:58:20.3950536Z building [mo]: targets for 0 po files that are out of date
2022-03-24T18:58:20.3968916Z building [html]: targets for 11 source files that are out of date
2022-03-24T18:58:20.3974478Z updating environment: [new config] 11 added, 0 changed, 0 removed
2022-03-24T18:58:20.3976429Z reading sources... [  9%] _examples/tutorial_1_wavebot
2022-03-24T18:58:20.9802350Z 
2022-03-24T18:58:20.9802956Z Notebook error:
2022-03-24T18:58:20.9803629Z AttributeError in _examples/tutorial_1_wavebot.ipynb:
2022-03-24T18:58:20.9804809Z module 'jinja2.utils' has no attribute 'escape'
2022-03-24T18:58:21.4410306Z make: *** [Makefile:22: html] Error 2
2022-03-24T18:58:21.4429650Z ##[error]Process completed with exit code 2.

CI/CD

Set this up please ;)

consolidate/remove plot_continuous

API documentation and class/method/function references should not use 'core' directly

Currently the API documentation refers to functions in core as e.g. wecopttool.core.WEC. Everything in core is imported into wecopttool in the __init__.py file, and the recommended use is e.g. wecopttool.WEC. This use should be reflected in the API documentation and any references to these functions in the rest of the documentation.
It shoud be possible through a combination of:
https://stackoverflow.com/questions/30856279/how-to-use-sphinx-automodule-and-exposed-functions-in-init
and specifying specific classes/functions to document rather than the entire wecopttool.

Move to Python v3.10

We are using functionality from Python 3.10 using from future import ... but cannot support python 3.10 yet because there are no VTK Python 3.10 wheels for PyPI yet. Once these are available we can either:

1. Move completely to Python>=3.10 . This would mean we don't need to use from future import ..., but would also mean users would need to have Python'>='3.10 (currently only 3.8 and 3.9 work).
2. Keep using from future import ... and change the Python requirement to >= 3.8.

Feature request: PTO with non-linear kinematics

Feature description.
Include a PTO class in the pto module that can handle non-linear kinematics.

Issue addressed.
Currently all PTO examples use linear kinematics.

Describe the solution you'd like
A new PTO class (or more) that can handle non-linear kinematics. Instead of passing a matrix the user would pass a function of the state.

Describe alternatives you've considered
Users can implement this themselves, but it would be good to include an example as is this is a common application.

Interest in leading this feature development?
Yes.

Additional information
The pto module will be restructured soon. For now implement based on current structure. Make sure this works with other PR's: linear PTO, PTO<=>WEC DoF.

BUG: WaveBot geom not using freeboard properly

in geom.py the function WaveBot.plot_cross_section plots the freeboard below the water line.

Improve documentation

The documentation should be improved. Highest priority right now:

• Explain the theory in more detail
• Explain how to choose scaling
• Explain why one might want to use more points for constraints than for the dynamics.

replicate case A from "Initial conceptual..."

Replicate case A from our previous paper (which used MATLAB, not Python). This compares three control structures/solutions:

• CC: complex conjugate*
• P: proportional damping*
• PS: pseudo-spectral (w. max PTO force of 2 kN)

*Note that in this paper CC was solved via analytic equations (see, e.g., sec 3.5 in Falnes) and P was solved by numerical optimization in the frequency domain. Both can be solved using the pseudo-spectral method. For the complex conjugate, we can simply remove the limit on PTO force. For the proportional damping, we should be able to use a constraint to prevent the controller from using reactive power.

Coe, R.G., Bacelli, G., Olson, S. et al. Initial conceptual demonstration of control co-design for WEC optimization. J. Ocean Eng. Mar. Energy 6, 441–449 (2020). https://doi.org/10.1007/s40722-020-00181-9

Feature request: Add linear force components to post-processing.

Feature description.
Calculate the different force components for the optimized solution (using the optimal state + the BEM coefficients) and include it in the xarray in the post-processing function. These forces include: excitation (split into diffraction + F-K?), radiation damping, hydrostatic, added mass.

Zero frequency component & intuitive buoyancy/gravity

In computation, we use an impedance where position is the state (flow) variable. This is called Gi in our parlance. Thus, we add a zero frequency component equal to the linear stiffness. This done when we create a block matrix.

The zero frequency component is not stored in wec.hydro, even though wec.hydro.Gi exists.

Our code runs, but is somewhat inconsistent and could probably be improved.

automatic type checking and types in docstrings

I was working on #78, specifically addressing #75... scipy.minimize.optimize can handle two types for the bounds argument: Bounds or sequence. This is tricky to handle, so I wrote code to handle only Bounds types for the bounds_wec and bounds_opt to wec.solve. That got me thinking about checking the types and raising a more helpful error than whatever happens randomly. I was about to do something like

if not isinstance(bounds_wec, Bounds):
    raise TypeError('...')

but then I read some more and it looks like there are some sleek solutions for doing this based entirely on the type hints that you put in the function signatures:

• typeguard: https://github.com/agronholm/typeguard
• beartype: https://github.com/beartype/beartype (from my reading, this looks like the better option)
• mypy: http://mypy-lang.org

If I understand correctly, we could use beartype by:

1. adding it to our dependencies

2. calling from beartype import beartype in each module

3. using the @beartype decorator on functions and methods

Any thoughts @cmichelenstrofer?

Bounds for decision variable

In certain cases, the user will want to set bounds for the decision variable elements. For example, when using a proportional damping controller the gain should always be negative (or positive, depending on the convention).

To accomplish this, we can take any additional keyword arguments for wec.solve and pass those along to scipy.optimize.minimize. This would also enable other more detail control over scipy.optimize.minimize in the future.

BUG: link check error from sphinx

Without a change to the source of our docs, we started getting a 403 error when running make clean linkcheck for https://doi.org/10.3390/en10040472

Broken link in RELEASING

Should point to main repo, not @cmichelenstrofer's fork

https://github.com/SNL-WaterPower/WecOptTool/blob/9eb7a90302c649a280a9ce6307fa6d44fd877e73/.github/RELEASING.md?plain=1#L27

Feature request: example with non-trivial buoyancy/gravity

Issue addressed.
capability to explicitly account for buoyancy, gravity, and pretension exists, but probably not too easy for external user to understand.

Describe the solution you'd like
add an example (or alter an existing one) to demonstrate this (e.g., a version of the WaveBot in which a pretension is applied to maintain equilibrium).

Describe alternatives you've considered
There is some language in the docs, but example would probably still be helpful.

Interest in leading this feature development?
Yes

BUG: initial guesses must be manually scaled

Describe the bug
When the user sets an initial guess for the decision variable (via x_opt_0 and/or x_wec_0) passed to WEC.solve and also uses the scale_x_wec, scale_x_opt, and scale_obj arguments, she must do something like x_opt_0 = my_guess * scale_x_opt

To Reproduce
Steps to reproduce the behavior:

1. Run wec.solve with scale_x_wec, scale_x_opt, and scale_obj arguments set
2. Not the value of x_opt3.
3. Rerun wec.solve, this time setting x_opt_0 with the value from step 2

Expected behavior
Very fast convergence

Observed behavior
Sometimes divergence (if scale_x_wec, scale_x_opt, and scale_obj are not close to unity)

System:

• OS: macOS
• Python version: 3.9.9
• WecOptTool version: 1.0.2

Additional information
@cmichelenstrofer suggest that this can be resolved by modifying the following line

https://github.com/SNL-WaterPower/WecOptTool/blob/27cf9b365af77511e14320de4f419fdfb489b024/wecopttool/core.py#L691

Feature request: Linear PTO

Feature description.
Include an example linear PTO in the pto module.

Issue addressed.
Currently example PTOs only look at mechanical power. A linear PTO will allow considering electric power.

Describe the solution you'd like
Create a new PTO class that can accept a PTO impedance and can return electric quantities (current, voltage, power).

Describe alternatives you've considered
The PTO module is more of an example repository, users can implement any PTO they like. But a linear PTO is a good baseline one to include as an example.

Interest in leading this feature development?
Yes.

Additional information
The PTO module will likely be restructured after this and other PTO's are implemented.

Feature request: cache deps. for CI/CD

https://docs.github.com/en/actions/advanced-guides/caching-dependencies-to-speed-up-workflows

BUG: freq. domain Dataset, power at wrong frequencies

While working on #38, we noticed that when plotting the power in the frequency domain, we see nonzero values where we don't expect them. This xr.Dataset is produced by taking fft's of the time domain results. The time domain results agree with analytical solutions, so we suspect something wrong with this post-processing workflow.

Wave directionality

Our current tensor structures are set up to handle multiple modes of motion, which may arise from a multi-modal device and/or multiple devices within an array. However, we do not currently have a means of handling directional waves. This should be simple enough however and relatively undisruptive to introduce, since the additional dimension needed to handle wave directionality will only show up when calculating the excitation force, after which this dimension is projected into the modes of motion.

Better scaling and/or initial guess

Do one or both of the following:

• Run w/o inequality constraints
• find the analytical solution (v = Fe/2/Ri)

From these, you can do one or both of the following:

• set x0
• set scaling

related to #4 #48

scaling

Since the forces and moments are often on the scale of 1e3 to 1e6 and velocities are on the scale 1e-1 to 1e0, the optimization problem needs some help with scaling... @cmichelenstrofer has already done this elsewhere, just need to apply here. Might be nice if you can make a note on this issue about the logic behind what you've done.