Warning: This document is for an old version of WecOptTool.

2. Example

This section explains and expands upon the example.m example file provided in the root directory of the WecOptTool source code. This example considers the DOE Reference Model 3 (RM3) device.

See the entire example file

% Copyright 2020 National Technology & Engineering Solutions of Sandia, 
% LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the 
% U.S. Government retains certain rights in this software.
%
% This file is part of WecOptTool.
% 
%     WecOptTool is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     WecOptTool is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with WecOptTool.  If not, see <https://www.gnu.org/licenses/>.

%% Suppress warnings
warning('off', 'WaveSpectra:NoWeighting')

%% Create an RM3 study object
study = WecOptTool.RM3Study();

%% Create Bretschnider spectrum from WAFO
%S = bretschneider([],[8,10],0);

%% Alternatively load a single example spectrum
% S = WecOptLib.tests.data.exampleSpectrum();

%% Or load an example with multiple sea-states (8 differing spectra)
S = WecOptLib.tests.data.example8Spectra();

%% Add spectra to study
study.addSpectra(S);

%% Add a controller to the study
cc = WecOptTool.control.ComplexConjugate();
study.addControl(cc);

%% Add geometry design variables (parametric)
x0 = [5, 7.5, 1.125, 42];
lb = [4.5, 7, 1.00, 41];
ub = [5.5, 8, 1.25, 43];

parametric = WecOptTool.geom.Parametric(x0, lb, ub);
study.addGeometry(parametric);

%% Add geometry design variables (scalar)
% x0 = 1.;
% lb = 0.5;
% ub = 2;
% 
% scalar = WecOptTool.geom.Scalar(x0, lb, ub);
% study.addGeometry(scalar);

%% Set up some optimisation options for fmincon
options = optimoptions('fmincon');
options.MaxFunctionEvaluations = 5; % set artificial low for fast running
options.UseParallel = true;

%% Run the study
WecOptTool.run(study, options);

%% Print the results
WecOptTool.result(study);

%% Plot the results
WecOptTool.plot(study);

%% Re-enable warnings
warning('on', 'WaveSpectra:NoWeighting')

The general concept of WecOptTool is illustrated in the diagram below. In the upper left-hand corner, an optimization algorithm controls the selection of a set of design variables. In this diagram, some geometric design variables, \(r_1, r_2, d_1, d_2\), are considered along with constraints on the power take-off (PTO) max force, \(F_{max}\), and max stroke \(\Delta x_{max}\) and an operational constraint, \(H_{s,max}\). The device defined by these design variables is passed to the grey evaluation block. Here, Nemoh is used to compute the linear wave-body interaction properties using the boundary element method (BEM). Next, using these properties and some set of sea states, one of three controllers (ProportionalDamping, ComplexConjugate, PseudoSpectral) are used to compute the resulting dynamics. From on these dynamics and some model for cost (e.g., based on the size dimensions and the capabilities of the PTO) can be combine to produce an objective function, which is returned to the optimization solver.

2.1. Create an RM3Study Object

The RM3Study class allows the user to configure a simulation to their specifications. Once instantiated, an RM3Study object can be modified using other classes (as described below), and once prepared is passed to the main functions of the toolbox.

warning('off', 'WaveSpectra:NoWeighting')

2.2. Define a Sea-State

WecOptTool can simulate single or multiple spectra sea states, where weightings can be provided to indicate the relative likelihood of each spectra. The following lines from example.m provide means of using the WAFO MATLAB toolbox or preset spectra from WecOptTool.

study = WecOptTool.RM3Study();

%% Create Bretschnider spectrum from WAFO
%S = bretschneider([],[8,10],0);

%% Alternatively load a single example spectrum
% S = WecOptLib.tests.data.exampleSpectrum();

Spectra are formatted following the convention of the WAFO MATLAB toolbox, but can be generated in via any means (e.g., from buoy measurements) as long as the structure includes the S.S, S.w, and S.phi fields.

S =

  struct with fields:

       S: [257×1 double]
       w: [257×1 double]
      tr: []
       h: Inf
    type: 'freq'
     phi: 0
    norm: 0
    note: 'Bretschneider, Hm0 = 4, Tp = 5'
    date: '25-Mar-2020 13:08:28'

In the active code above from example.m, there are eight spectra loaded into a struct array. These can be plotted using standard MATLAB commands.

figure
hold on
grid on
arrayfun(@(x) plot(x.w,x.S,'DisplayName',x.note), S)
legend()
xlim([0,3])
xlabel('Freq. [rad/s]')
ylabel('Spect. density [m^2 rad/s]')

The desired spectrum or spectra can then be added to the study object.

S = WecOptLib.tests.data.example8Spectra();

2.3. Add a controller to the study

WecOptTool allows for three types of controllers:

• ProportionalDamping: Resistive damping (i.e., a proportional feedback on velocity) (see, e.g., [Falnes]). Here, the power take-off (PTO) force is set as

  \[F_u(\omega) = -B_{PTO}(\omega)u(\omega)\]

  where \(B_{PTO}\) is a constant chosen to maximize absorbed power and \(u(\omega)\) is the velocity.
  
  

• ComplexConjugate: Optimal power absorption through impedance matching (see, e.g., [Falnes]). The intrinsic impedance is given by

  \[Z_i(\omega) = B(\omega) + i \left( \omega{}(m + A(\omega)) - \frac{K_{HS}}{\omega}\right) ,\]

  where \(\omega\) is the radial frequency, \(B(\omega)\) is the radiation damping, \(m\) is the rigid body mass, \(A(\omega)\) is the added mass, and \(K_{HS}\) is the hydrostatic stiffness. Optimal power transfer occurs when the PTO force, \(F_u\) is set such that

  \[F_u(\omega) = -Z_i^*(\omega)u(\omega) ,\]

  where \(u(\omega)\) is the velocity.
  
  

• PseudoSpectral: Constrained optimal power absorption [Bacelli]. This is a numerical optimal control algorithm capable of dealing with both constraints and nonlinear dynamics. This approach is based on pseudo spectral optimal control.

The controllers are defined as classes in the control sub-package.

study.addSpectra(S);

%% Add a controller to the study

2.4. Define design variables

As shown in the diagram below, for RM3 study considered in example.m the design variables are the radius of the surface float, r1, the radius of the heave plate, r2, the draft of the surface float, d1, and the depth of the heave plate, d2, such that x = [r1, r2, d1, d2]. The optimization algorithm will attempt to find the values of x that minimize the objective function.

Note

Objective function: The built-in objective function of example.m is set to maximize absorbed power. This function can be altered better approximate a more meaningful objective (e.g., levelized cost of energy).

The initial values,``x0``, lower bounds, lb, and upper bounds, ub of the design variables can be set as follows.

study.addControl(cc);

%% Add geometry design variables (parametric)
x0 = [5, 7.5, 1.125, 42];
lb = [4.5, 7, 1.00, 41];
ub = [5.5, 8, 1.25, 43];

Alternatively, a simpler study with a single scalar design variable can be employed. In this case, instead of scaling various dimensions of the device individually, the entire device is scaled based on a single design variable.

study.addGeometry(parametric);

%% Add geometry design variables (scalar)
% x0 = 1.;
% lb = 0.5;
% ub = 2;
% 

The options for design variables are defined as classes in the geom sub-package.

2.5. Set optimization solver and options

MATLAB’s fmincon optimization solver is used in example.m.

Note

The MaxFunctionEvaluations is set to 5 in example.m to permit relatively quick runs, but can be increased to allow for a potentially better solution (with the other options left as-is, this should require 150 function evaluations).

% study.addGeometry(scalar);

%% Set up some optimisation options for fmincon
options = optimoptions('fmincon');

2.6. Run the study and view results

The study can be run() and reviewed (with result() and plot()) as follows:

options.UseParallel = true;

%% Run the study
WecOptTool.run(study, options);

%% Print the results
WecOptTool.result(study);

From result(), we can see that the study has produced the following design.

Optimal solution is:
r1: 5 [m]
r2: 7.5 [m]
d1: 1.125 [m]
d2: 42 [m]
Optimal function value is: 2202421.2193 [W]

From plot(), we can review the spectral distribution of energy absorbed by the resulting design for each of the eight sea states.

[Falnes]

(1, 2) Falnes, Johannes. Ocean waves and oscillating systems: linear interactions including wave-energy extraction. Cambridge University Press, 2002.

[Bacelli]

Bacelli, Giorgio, and John V. Ringwood. “Numerical optimal control of wave energy converters.” IEEE Transactions on Sustainable Energy 6.2 (2014): 294-302.