From ac9a9f4be8322c56205a34c49e6e5ac448dbdbd2 Mon Sep 17 00:00:00 2001 From: hivanov-nrel Date: Tue, 23 Sep 2025 12:13:13 -0600 Subject: [PATCH 01/12] bug fix for phase angle --- mhkit/loads/extreme/mler_coefficients.m | 2 +- mhkit/tests/Loads_TestExtreme.m | 8 +------- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 165b83358..e659e73ed 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -52,7 +52,7 @@ % calculate coefficient_a from Quon2016 Eqn.8 coeff_a_rn = abs(RAO) .* sqrt(2*wave_spectrum*dw) .* ((m2 - freq*m1) + wBar*(freq*m0 - m1)) ./ (m0*m2 - m1^2); % phase delay should be positive number -phase = -unwrap(angle(RAO)); +phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign phase(coeff_a_rn < 0) = phase(coeff_a_rn < 0) - pi; coeff_a_rn(coeff_a_rn < 0) = coeff_a_rn(coeff_a_rn < 0) * -1; diff --git a/mhkit/tests/Loads_TestExtreme.m b/mhkit/tests/Loads_TestExtreme.m index f6ef72099..dd16a1be8 100644 --- a/mhkit/tests/Loads_TestExtreme.m +++ b/mhkit/tests/Loads_TestExtreme.m @@ -3,17 +3,11 @@ methods (Test) function test_mler_coefficients(testCase) - % `mler_coefficients` Line 54, 55 - % phase delay should be positive number - % phase = -unwrap(angle(RAO)); - % Per @simmsa, this does not perform this conversion correctly for positive RAO values - assumeFail(testCase, "Per @simmsa, ask @hivanov about setting RAO to negative values per the above comment"); - % create inputs and load validation data fpath = '../../examples/data/loads/mler.csv'; validation = readtable(fpath); wave_freq = linspace(0,1,500); - js = jonswap_spectrum(wave_freq,15.1,9); + js = pierson_moskowitz_spectrum(wave_freq,15.1,9); response_desired = 1; RAO = validation.RAO; RAO = RAO'; From 542301f523623bd0a3fc7a10e665f286573346e1 Mon Sep 17 00:00:00 2001 From: hivanov-nrel Date: Tue, 23 Sep 2025 12:38:40 -0600 Subject: [PATCH 02/12] explicit elmntwise multiplication --- mhkit/loads/extreme/mler_coefficients.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index e659e73ed..01c52b157 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -50,7 +50,7 @@ wBar = m1/m0; % calculate coefficient_a from Quon2016 Eqn.8 -coeff_a_rn = abs(RAO) .* sqrt(2*wave_spectrum*dw) .* ((m2 - freq*m1) + wBar*(freq*m0 - m1)) ./ (m0*m2 - m1^2); +coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); % phase delay should be positive number phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign From 7ec6ead9ae8929e9ded0ab51afcbef4380a4b6e9 Mon Sep 17 00:00:00 2001 From: hivanov-nrel Date: Tue, 23 Sep 2025 12:53:37 -0600 Subject: [PATCH 03/12] multiply fix v2 --- mhkit/loads/extreme/mler_coefficients.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 01c52b157..3674e943e 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -50,7 +50,7 @@ wBar = m1/m0; % calculate coefficient_a from Quon2016 Eqn.8 -coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); +coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0.*m2 - m1.^2); % phase delay should be positive number phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign From 69b8872fecce8f563a42463ceb5642be7309a6b4 Mon Sep 17 00:00:00 2001 From: hivanov-nrel Date: Tue, 23 Sep 2025 15:09:16 -0600 Subject: [PATCH 04/12] multiply fix v3 --- mhkit/loads/extreme/mler_coefficients.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 3674e943e..01c52b157 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -50,7 +50,7 @@ wBar = m1/m0; % calculate coefficient_a from Quon2016 Eqn.8 -coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0.*m2 - m1.^2); +coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); % phase delay should be positive number phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign From c4a990f63319004c061c1e78b9a128f8cc686716 Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 16:23:08 -0600 Subject: [PATCH 05/12] Tests: Explicitly print types in mler_coefficients --- mhkit/loads/extreme/mler_coefficients.m | 61 ++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 01c52b157..808b0c7fc 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -49,8 +49,67 @@ m2 = frequency_moment(R, 2); wBar = m1/m0; +% DEBUG: Print all dimensions before calculation +fprintf('DEBUG - Variable dimensions in mler_coefficients:\n'); +fprintf(' RAO: %s\n', mat2str(size(RAO))); +fprintf(' wave_spectrum: %s\n', mat2str(size(wave_spectrum))); +fprintf(' freq: %s\n', mat2str(size(freq))); +fprintf(' dw: %s (scalar: %g)\n', mat2str(size(dw)), dw); +fprintf(' m0: %s (value: %g)\n', mat2str(size(m0)), m0); +fprintf(' m1: %s (value: %g)\n', mat2str(size(m1)), m1); +fprintf(' m2: %s (value: %g)\n', mat2str(size(m2)), m2); +fprintf(' wBar: %s (value: %g)\n', mat2str(size(wBar)), wBar); + % calculate coefficient_a from Quon2016 Eqn.8 -coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); +% Explicit scalar handling for MATLAB R2022b/R2023b compatibility +% Debug output before computation +fprintf(' BEFORE explicit scalar conversion:\n'); +fprintf(' m0 original: size=%s, class=%s, value=%g\n', mat2str(size(m0)), class(m0), m0); +fprintf(' m1 original: size=%s, class=%s, value=%g\n', mat2str(size(m1)), class(m1), m1); +fprintf(' m2 original: size=%s, class=%s, value=%g\n', mat2str(size(m2)), class(m2), m2); +fprintf(' wBar original: size=%s, class=%s, value=%g\n', mat2str(size(wBar)), class(wBar), wBar); + +% Ensure all scalars are explicitly converted to same size as vectors +m0_scalar = double(m0(1)); % Ensure scalar +m1_scalar = double(m1(1)); % Ensure scalar +m2_scalar = double(m2(1)); % Ensure scalar +wBar_scalar = double(wBar(1)); % Ensure scalar + +fprintf(' AFTER explicit scalar conversion:\n'); +fprintf(' m0_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m0_scalar)), class(m0_scalar), m0_scalar); +fprintf(' m1_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m1_scalar)), class(m1_scalar), m1_scalar); +fprintf(' m2_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m2_scalar)), class(m2_scalar), m2_scalar); +fprintf(' wBar_scalar: size=%s, class=%s, value=%g\n', mat2str(size(wBar_scalar)), class(wBar_scalar), wBar_scalar); + +% Calculate each term explicitly with size monitoring +term1 = abs(RAO); % [1 x N] +fprintf(' term1 = abs(RAO): %s\n', mat2str(size(term1))); + +term2 = sqrt(2 .* dw .* wave_spectrum); % [1 x N] +fprintf(' term2 = sqrt(2.*dw.*wave_spectrum): %s\n', mat2str(size(term2))); + +term3a = m2_scalar - (freq .* m1_scalar); % [1 x N] +fprintf(' term3a = m2_scalar - (freq.*m1_scalar): %s\n', mat2str(size(term3a))); + +term3b = (freq .* m0_scalar) - m1_scalar; % [1 x N] +fprintf(' term3b = (freq.*m0_scalar) - m1_scalar: %s\n', mat2str(size(term3b))); + +term3c = wBar_scalar .* term3b; % [1 x N] +fprintf(' term3c = wBar_scalar .* term3b: %s\n', mat2str(size(term3c))); + +term3 = term3a + term3c; % [1 x N] +fprintf(' term3 = term3a + term3c: %s\n', mat2str(size(term3))); + +% Final calculation with explicit element-wise operations +numerator = term1 .* term2 .* term3; % [1 x N] +fprintf(' numerator = term1 .* term2 .* term3: %s\n', mat2str(size(numerator))); + +denominator_scalar = (m0_scalar * m2_scalar) - (m1_scalar ^ 2); % scalar +fprintf(' denominator_scalar: size=%s, value=%g\n', mat2str(size(denominator_scalar)), denominator_scalar); + +coeff_a_rn = numerator ./ denominator_scalar; % [1 x N] +fprintf(' coeff_a_rn final: %s\n', mat2str(size(coeff_a_rn))); + % phase delay should be positive number phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign From 4b24e0d47a812a83aa38a6a354e12e9225e096a8 Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 16:33:52 -0600 Subject: [PATCH 06/12] Revert "Tests: Explicitly print types in mler_coefficients" This reverts commit c4a990f63319004c061c1e78b9a128f8cc686716. --- mhkit/loads/extreme/mler_coefficients.m | 61 +------------------------ 1 file changed, 1 insertion(+), 60 deletions(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 808b0c7fc..01c52b157 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -49,67 +49,8 @@ m2 = frequency_moment(R, 2); wBar = m1/m0; -% DEBUG: Print all dimensions before calculation -fprintf('DEBUG - Variable dimensions in mler_coefficients:\n'); -fprintf(' RAO: %s\n', mat2str(size(RAO))); -fprintf(' wave_spectrum: %s\n', mat2str(size(wave_spectrum))); -fprintf(' freq: %s\n', mat2str(size(freq))); -fprintf(' dw: %s (scalar: %g)\n', mat2str(size(dw)), dw); -fprintf(' m0: %s (value: %g)\n', mat2str(size(m0)), m0); -fprintf(' m1: %s (value: %g)\n', mat2str(size(m1)), m1); -fprintf(' m2: %s (value: %g)\n', mat2str(size(m2)), m2); -fprintf(' wBar: %s (value: %g)\n', mat2str(size(wBar)), wBar); - % calculate coefficient_a from Quon2016 Eqn.8 -% Explicit scalar handling for MATLAB R2022b/R2023b compatibility -% Debug output before computation -fprintf(' BEFORE explicit scalar conversion:\n'); -fprintf(' m0 original: size=%s, class=%s, value=%g\n', mat2str(size(m0)), class(m0), m0); -fprintf(' m1 original: size=%s, class=%s, value=%g\n', mat2str(size(m1)), class(m1), m1); -fprintf(' m2 original: size=%s, class=%s, value=%g\n', mat2str(size(m2)), class(m2), m2); -fprintf(' wBar original: size=%s, class=%s, value=%g\n', mat2str(size(wBar)), class(wBar), wBar); - -% Ensure all scalars are explicitly converted to same size as vectors -m0_scalar = double(m0(1)); % Ensure scalar -m1_scalar = double(m1(1)); % Ensure scalar -m2_scalar = double(m2(1)); % Ensure scalar -wBar_scalar = double(wBar(1)); % Ensure scalar - -fprintf(' AFTER explicit scalar conversion:\n'); -fprintf(' m0_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m0_scalar)), class(m0_scalar), m0_scalar); -fprintf(' m1_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m1_scalar)), class(m1_scalar), m1_scalar); -fprintf(' m2_scalar: size=%s, class=%s, value=%g\n', mat2str(size(m2_scalar)), class(m2_scalar), m2_scalar); -fprintf(' wBar_scalar: size=%s, class=%s, value=%g\n', mat2str(size(wBar_scalar)), class(wBar_scalar), wBar_scalar); - -% Calculate each term explicitly with size monitoring -term1 = abs(RAO); % [1 x N] -fprintf(' term1 = abs(RAO): %s\n', mat2str(size(term1))); - -term2 = sqrt(2 .* dw .* wave_spectrum); % [1 x N] -fprintf(' term2 = sqrt(2.*dw.*wave_spectrum): %s\n', mat2str(size(term2))); - -term3a = m2_scalar - (freq .* m1_scalar); % [1 x N] -fprintf(' term3a = m2_scalar - (freq.*m1_scalar): %s\n', mat2str(size(term3a))); - -term3b = (freq .* m0_scalar) - m1_scalar; % [1 x N] -fprintf(' term3b = (freq.*m0_scalar) - m1_scalar: %s\n', mat2str(size(term3b))); - -term3c = wBar_scalar .* term3b; % [1 x N] -fprintf(' term3c = wBar_scalar .* term3b: %s\n', mat2str(size(term3c))); - -term3 = term3a + term3c; % [1 x N] -fprintf(' term3 = term3a + term3c: %s\n', mat2str(size(term3))); - -% Final calculation with explicit element-wise operations -numerator = term1 .* term2 .* term3; % [1 x N] -fprintf(' numerator = term1 .* term2 .* term3: %s\n', mat2str(size(numerator))); - -denominator_scalar = (m0_scalar * m2_scalar) - (m1_scalar ^ 2); % scalar -fprintf(' denominator_scalar: size=%s, value=%g\n', mat2str(size(denominator_scalar)), denominator_scalar); - -coeff_a_rn = numerator ./ denominator_scalar; % [1 x N] -fprintf(' coeff_a_rn final: %s\n', mat2str(size(coeff_a_rn))); - +coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); % phase delay should be positive number phase = unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign From 4f4e4e022f32c60ad6690f264c33cef0da8e78a3 Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 17:18:15 -0600 Subject: [PATCH 07/12] Loads: Specify input dimensions in mler_coefficients docstring --- mhkit/loads/extreme/mler_coefficients.m | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index 01c52b157..cf73c8d25 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -7,18 +7,20 @@ % % Parameters % ---------- -% RAO : array +% RAO : array (N x 1) % Response amplitude operator [-] % wave_spectrum: struct -% Struct with wave spectral density [m^2/Hz] and frequency [Hz] -% response_desired: int or float -% Latitude longitude pairs at which to extract data. +% wave_spectrum.spectrum - Spectral density [m^2/Hz] (N x 1) +% wave_spectrum.frequency - Frequency [Hz] (N x 1) +% response_desired: scalar +% Desired response amplitude % % Returns % ------- % mler : struct -% containing conditioned wave spectral amplitude -% coefficient [m^2-s], and phase [rad] indexed by frequency [Hz]. +% mler.conditioned_spectrum - Conditioned wave spectral amplitude [m^2-s] (N x 1) +% mler.phase - Phase [rad] (N x 1) +% mler.frequency - Frequency [Hz] (N x 1) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% From 55513bddfffb2ea2ef3d58ae8e17ca21a767edbe Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 17:19:09 -0600 Subject: [PATCH 08/12] Loads: Be explicit about how column vs row vectors are handled/validated --- mhkit/loads/extreme/mler_coefficients.m | 29 +++++++++++++++---------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/mhkit/loads/extreme/mler_coefficients.m b/mhkit/loads/extreme/mler_coefficients.m index cf73c8d25..7fed12248 100644 --- a/mhkit/loads/extreme/mler_coefficients.m +++ b/mhkit/loads/extreme/mler_coefficients.m @@ -34,16 +34,21 @@ error('ERROR: response_desired must be an int or double') end +% validate all inputs have same length +N = length(RAO); +if length(wave_spectrum.frequency) ~= N || length(wave_spectrum.spectrum) ~= N + error('MHKiT:loads:mler_coefficients: RAO, frequency, and spectrum must have same length'); +end + % convert from Hz to rad/s -freq = wave_spectrum.frequency * (2*pi); -freq_hz = wave_spectrum.frequency; -wave_spectrum = wave_spectrum.spectrum / (2*pi); -dw = (2*pi - 0) / (length(freq)-1); +freq_rad = wave_spectrum.frequency * (2*pi); +wave_spectrum_rad = wave_spectrum.spectrum / (2*pi); +dw = (2*pi - 0) / (N-1); % response spectrum -R.spectrum = abs(RAO).^2 .* (2*wave_spectrum); +R.spectrum = abs(RAO).^2 .* (2*wave_spectrum_rad); R.type = 'response'; -R.frequency = freq; +R.frequency = freq_rad; % spectral moment calculations m0 = frequency_moment(R, 0); @@ -52,16 +57,16 @@ wBar = m1/m0; % calculate coefficient_a from Quon2016 Eqn.8 -coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2); +coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum_rad) .* ((m2 - freq_rad.*m1) + wBar.*(freq_rad.*m0 - m1)) ./ (m0*m2 - m1^2); % phase delay should be positive number -phase = unwrap(angle(RAO)); +phase = -unwrap(angle(RAO)); % for negative values of Amp, add pi phase shift, flip sign phase(coeff_a_rn < 0) = phase(coeff_a_rn < 0) - pi; coeff_a_rn(coeff_a_rn < 0) = coeff_a_rn(coeff_a_rn < 0) * -1; % calculate conditioned spectrum [m^2-s/rad] -S = wave_spectrum .* coeff_a_rn.^2 .* response_desired^2; -S(isnan(S)) = 0; % replace nans with zero +conditioned_spectrum = wave_spectrum_rad .* coeff_a_rn.^2 .* response_desired^2; +conditioned_spectrum(isnan(conditioned_spectrum)) = 0; % replace nans with zero % if the response amplitude we ask for is negative, we will add % a pi phase shift to the phase information. This is because % the sign of response_desired is lost in the squaring above. @@ -73,8 +78,8 @@ end % outputs -mler.conditioned_spectrum = S; +mler.conditioned_spectrum = conditioned_spectrum; mler.phase = phase; -mler.frequency = freq_hz; +mler.frequency = wave_spectrum.frequency; end From 0aaf524f87bd11e584b6f094404b8996eeecc4e0 Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 17:20:20 -0600 Subject: [PATCH 09/12] Tests: Remove transposing from test_mler_coefficients --- mhkit/tests/Loads_TestExtreme.m | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/mhkit/tests/Loads_TestExtreme.m b/mhkit/tests/Loads_TestExtreme.m index dd16a1be8..603fa6701 100644 --- a/mhkit/tests/Loads_TestExtreme.m +++ b/mhkit/tests/Loads_TestExtreme.m @@ -10,12 +10,11 @@ function test_mler_coefficients(testCase) js = pierson_moskowitz_spectrum(wave_freq,15.1,9); response_desired = 1; RAO = validation.RAO; - RAO = RAO'; % execute function mler = mler_coefficients(RAO, js, response_desired); % assertions - assertEqual(testCase, mler.conditioned_spectrum, validation.Res_Spec', 'RelTol',0.005) - assertEqual(testCase, mler.phase, validation.phase', 'RelTol',0.001) + assertEqual(testCase, mler.conditioned_spectrum, validation.Res_Spec, 'RelTol',0.005) + assertEqual(testCase, mler.phase, validation.phase, 'RelTol',0.001) end function test_mler_simulation(testCase) From d7a227c3f9350c242d4e0feb974d3042c3e6b80a Mon Sep 17 00:00:00 2001 From: "Simms, Andrew" Date: Tue, 23 Sep 2025 17:50:55 -0600 Subject: [PATCH 10/12] Examples: MLER, remove unnecessary transpose --- examples/extreme_response_MLER_example.html | 11 +++-------- examples/extreme_response_MLER_example.m | 7 +------ examples/extreme_response_MLER_example.mlx | Bin 203056 -> 103064 bytes 3 files changed, 4 insertions(+), 14 deletions(-) diff --git a/examples/extreme_response_MLER_example.html b/examples/extreme_response_MLER_example.html index e35738089..175a5eda5 100644 --- a/examples/extreme_response_MLER_example.html +++ b/examples/extreme_response_MLER_example.html @@ -40,7 +40,7 @@ .inlineElement .textElement {} .embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement { min-height: 16px;} .rightPaneElement .textElement { padding-top: 2px; padding-left: 9px;} -.S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

MLER example

Extreme conditions modeling consists of identifying the expected extreme (e.g. 100-year) response of some quantity of interest, such as WEC motions or mooring loads. Three different methods of estimating extreme conditions were adapted from WDRT: full sea state approach, contour approach, and MLER design wave. This livescript presents the MLER approach.
This example notebook shows users how to utilize the most likely extreme response (MLER) method. This method is an alternative to exhaustive Monte Carlo or long-term simulations for finding and evaluating wave energy converter response events at extreme loads. To accomplish this, a given RAO is combined with a wave spectrum corresponding to certain extreme sea states. MLER then generates the resulting wave profile that would cause the largest response for that degree of freedom. See full explanation and derivation of the MLER method in
E. Quon, A. Platt, Y.-H. Yu, and M. Lawson, “Application of the Most Likely Extreme Response Method for Wave Energy Converters,” in Volume 6: Ocean Space Utilization; Ocean Renewable Energy, Busan, South Korea, Jun. 2016, p. V006T09A022. doi: 10.1115/OMAE2016-54751.
In this example, a simple ellipsoid shaped WEC device was modeled in WEC-Sim. We will focus on analyzing the heave response of this device. The code below simply imports RAO data as it is needed for one of the inputs.
wave_freq = linspace(0, 1, 500);
mfile = readtable('data/loads/mler.csv');
RAO = mfile.RAO;
Next, we need to generate a wave environment that corresponds to a chosen extreme sea state. The associated parameters are selected in different ways. In this case, public wave data was used to come up with estimated 100-year sea state contour. The sea state parameters were selected somewhere along the contour.
Hs = 9.0; % significant wave height
Tp = 15.1; % time period of waves
js = jonswap_spectrum(wave_freq, Tp, Hs);
 
figure('Position', [100, 100, 1600, 600]);
plot(js.frequency, js.spectrum); xlabel('frequency [Hz]'); ylabel('response [m^2/Hz]')
Now that we have both the RAO and the spectrum of the wave environment, we can calculate the MLER conditioned wave spectrum and phase. In this case, we would like to find the wave profile that will generate a heave response of 1 meter for our WEC device.
mler_data = mler_coefficients(RAO', js, 1);
plot(mler_data.frequency, mler_data.conditioned_spectrum); xlabel('Frequency [Hz]'); ylabel('Conditoned wave spectrum [m^2-s]')
plot(mler_data.frequency, mler_data.phase); xlabel('Frequency [Hz]'); ylabel('Phase [rad]')
From here, we can choose to export these coefficients to feed into other high fidelity software. However, if we wanted to get a specific height of the incoming wave, we can renormalize the wave amplitude. To do this, several inputs need to be generated in addition to the existing coefficients.
The first is a struct containing information pertinent to creating time series. This struct can be easily generated using mler_simulation. In this example, the input dictionary contains the default values of this function, but the format is shown in case a user wants to adjust the parameters for their specific use case. The second input is the wave number, which was obtained using the wave_number function from the resource module. Finally, we decide to renormalize the wave amplitude to a desired peak height (peak to MSL). Combining all of these inputs into mler_wave_amp_normalize gives us our new normalized mler wave spectrum.
sim.startTime = -150.0; % .s Starting time
sim.endTime = 150.0; % .s Ending time
sim.dT = 1.0; % .s Time-step size
sim.T0 = 0.0; % .s Time of maximum event
 
sim.startX = -300.0; % .m Start of simulation space
sim.endX = 300.0; % .m End of simulation space
sim.dX = 1.0; % .m Horizontal spacing
sim.X0 = 0.0; % .m Position of maximum event
 
sim = mler_simulation(sim)
sim = struct with fields:
startTime: -150 +.S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

MLER example

Extreme conditions modeling consists of identifying the expected extreme (e.g. 100-year) response of some quantity of interest, such as WEC motions or mooring loads. Three different methods of estimating extreme conditions were adapted from WDRT: full sea state approach, contour approach, and MLER design wave. This livescript presents the MLER approach.
This example notebook shows users how to utilize the most likely extreme response (MLER) method. This method is an alternative to exhaustive Monte Carlo or long-term simulations for finding and evaluating wave energy converter response events at extreme loads. To accomplish this, a given RAO is combined with a wave spectrum corresponding to certain extreme sea states. MLER then generates the resulting wave profile that would cause the largest response for that degree of freedom. See full explanation and derivation of the MLER method in
E. Quon, A. Platt, Y.-H. Yu, and M. Lawson, “Application of the Most Likely Extreme Response Method for Wave Energy Converters,” in Volume 6: Ocean Space Utilization; Ocean Renewable Energy, Busan, South Korea, Jun. 2016, p. V006T09A022. doi: 10.1115/OMAE2016-54751.
In this example, a simple ellipsoid shaped WEC device was modeled in WEC-Sim. We will focus on analyzing the heave response of this device. The code below simply imports RAO data as it is needed for one of the inputs.
wave_freq = linspace(0, 1, 500);
mfile = readtable('data/loads/mler.csv');
RAO = mfile.RAO;
Next, we need to generate a wave environment that corresponds to a chosen extreme sea state. The associated parameters are selected in different ways. In this case, public wave data was used to come up with estimated 100-year sea state contour. The sea state parameters were selected somewhere along the contour.
Hs = 9.0; % significant wave height
Tp = 15.1; % time period of waves
js = jonswap_spectrum(wave_freq, Tp, Hs);
 
figure('Position', [100, 100, 1600, 600]);
plot(js.frequency, js.spectrum); xlabel('frequency [Hz]'); ylabel('response [m^2/Hz]')
Now that we have both the RAO and the spectrum of the wave environment, we can calculate the MLER conditioned wave spectrum and phase. In this case, we would like to find the wave profile that will generate a heave response of 1 meter for our WEC device.
mler_data = mler_coefficients(RAO, js, 1);
plot(mler_data.frequency, mler_data.conditioned_spectrum); xlabel('Frequency [Hz]'); ylabel('Conditoned wave spectrum [m^2-s]')
plot(mler_data.frequency, mler_data.phase); xlabel('Frequency [Hz]'); ylabel('Phase [rad]')
From here, we can choose to export these coefficients to feed into other high fidelity software. However, if we wanted to get a specific height of the incoming wave, we can renormalize the wave amplitude. To do this, several inputs need to be generated in addition to the existing coefficients.
The first is a struct containing information pertinent to creating time series. This struct can be easily generated using mler_simulation. In this example, the input dictionary contains the default values of this function, but the format is shown in case a user wants to adjust the parameters for their specific use case. The second input is the wave number, which was obtained using the wave_number function from the resource module. Finally, we decide to renormalize the wave amplitude to a desired peak height (peak to MSL). Combining all of these inputs into mler_wave_amp_normalize gives us our new normalized mler wave spectrum.
sim.startTime = -150.0; % .s Starting time
sim.endTime = 150.0; % .s Ending time
sim.dT = 1.0; % .s Time-step size
sim.T0 = 0.0; % .s Time of maximum event
 
sim.startX = -300.0; % .m Start of simulation space
sim.endX = 300.0; % .m End of simulation space
sim.dX = 1.0; % .m Horizontal spacing
sim.X0 = 0.0; % .m Position of maximum event
 
sim = mler_simulation(sim)
sim = struct with fields:
startTime: -150 endTime: 150 dT: 1 T0: 0 @@ -52,7 +52,7 @@ T: [-150 -149 -148 -147 -146 -145 -144 -143 -142 -141 -140 -139 -138 -137 -136 -135 -134 -133 -132 -131 -130 -129 -128 -127 -126 -125 -124 -123 -122 -121 -120 -119 -118 -117 -116 -115 -114 -113 -112 -111 -110 -109 -108 -107 … ] (1×301 double) maxIX: 601 X: [-300 -299 -298 -297 -296 -295 -294 -293 -292 -291 -290 -289 -288 -287 -286 -285 -284 -283 -282 -281 -280 -279 -278 -277 -276 -275 -274 -273 -272 -271 -270 -269 -268 -267 -266 -265 -264 -263 -262 -261 -260 -259 -258 -257 … ] (1×601 double) -
 
k = wave_number(wave_freq, 70);
k.values = fillmissing(k.values,'constant',0);
 
peakHeightDesired = Hs/2 * 1.9;
 
% Transpose MLER output
mler_data.conditioned_spectrum = mler_data.conditioned_spectrum';
mler_data.phase = mler_data.phase';
mler_data.frequency = mler_data.frequency';
 
mler_norm = mler_wave_amp_normalize(peakHeightDesired, mler_data, sim, k.values);
As a final step, a user might need to convert the MLER coefficients into a time series for input into higher fidelity software. We can do this by using the mler_export_time_series function. The result is a struct containing the wave height [m] and the linear response [*] indexed by time.
mler_ts = mler_export_time_series(RAO, mler_norm, sim, k.values);
 
figure('Position', [100, 100, 1600, 600]);
plot(mler_ts.time, mler_ts.linear_response, mler_ts.time, mler_ts.wave_height);
xlabel('Time (s)'); ylabel('[*] / [m]'); grid on; legend('Linear Response', 'Wave Height')
+
 
k = wave_number(wave_freq, 70);
k.values = fillmissing(k.values,'constant',0);
 
peakHeightDesired = Hs/2 * 1.9;
 
mler_norm = mler_wave_amp_normalize(peakHeightDesired, mler_data, sim, k.values);
As a final step, a user might need to convert the MLER coefficients into a time series for input into higher fidelity software. We can do this by using the mler_export_time_series function. The result is a struct containing the wave height [m] and the linear response [*] indexed by time.
mler_ts = mler_export_time_series(RAO, mler_norm, sim, k.values);
 
figure('Position', [100, 100, 1600, 600]);
plot(mler_ts.time, mler_ts.linear_response, mler_ts.time, mler_ts.wave_height);
xlabel('Time (s)'); ylabel('[*] / [m]'); grid on; legend('Linear Response', 'Wave Height')