<jupyter_start><jupyter_text>Extreme Conditions Modeling - Full Sea State ApproachExtreme 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](https://github.com/WEC-Sim/WDRT): full sea state approach, contour approach, and MLER design wave. This noteboook presents the full sea state approach. The full sea state approach consists of the following steps: 1. Take $N$ samples to represent the sea state. Each sample represents a small area of the sea state and consists of a representative $(H_{s}, T_{e})_i$ pair and a weight $W_i$ associated with the probability of that sea state area. 2. For each sample $(H_{s}, T_{e})_i$ calculate the short-term (e.g. 3-hours) extreme for the quantity of interest (e.g. WEC motions or mooring tension).3. Integrate over the entire sea state to obtain the long-term extreme. This is a sum of the products of the weight of each sea state times the short-term extreme. See more details and equations in> [1] Coe, <NAME>., <NAME>, <NAME>, and <NAME>. 2018. “Full Long-Term Design Response Analysis of a Wave Energy Converter.” Renewable Energy 116: 356–66.**NOTE:** Prior to running this example it is recommended to become familiar with `environmental_contours_example.ipynb` and `short_term_extremes_example.ipynb` since some code blocks are adapted from those examples and used here without the additional description. We start by importing the relevant modules, including `waves.contours` submodule which includes the sampling function, and `loads.extreme` which inlcudes the short-term extreme and full sea state integration functions.<jupyter_code>from mhkit.wave import resource, contours, graphics
from mhkit.loads import extreme
import matplotlib.pyplot as plt
from mhkit.wave.io import ndbc
import pandas as pd
import numpy as np<jupyter_output><empty_output><jupyter_text>Obtain and Process NDBC Buouy DataThe first step will be obtaining the environmental data and creating the contours.See `environmental_contours_example.ipynb` for more details and explanations of how this is being done in the following code block.<jupyter_code>parameter = 'swden'
buoy_number = '46022'
ndbc_available_data = ndbc.available_data(parameter, buoy_number)

years_of_interest = ndbc_available_data[ndbc_available_data.year < 2013]

filenames = years_of_interest['filename']
ndbc_requested_data = ndbc.request_data(parameter, filenames)

ndbc_data = {}
for year in ndbc_requested_data:
    year_data = ndbc_requested_data[year]
    ndbc_data[year] = ndbc.to_datetime_index(parameter, year_data)

Hm0_list = []
Te_list = []

# Iterate over each year and save the result in the initalized dictionary
for year in ndbc_data:
    year_data = ndbc_data[year]
    Hm0_list.append(resource.significant_wave_height(year_data.T))
    Te_list.append(resource.energy_period(year_data.T))

# Concatenate list of Series into a single DataFrame
Te = pd.concat(Te_list, axis=0)
Hm0 = pd.concat(Hm0_list, axis=0)
Hm0_Te = pd.concat([Hm0, Te], axis=1)

# Drop any NaNs created from the calculation of Hm0 or Te
Hm0_Te.dropna(inplace=True)
# Sort the DateTime index
Hm0_Te.sort_index(inplace=True)

Hm0_Te_clean = Hm0_Te[Hm0_Te.Hm0 < 20]

Hm0 = Hm0_Te_clean.Hm0.values
Te = Hm0_Te_clean.Te.values

dt = (Hm0_Te_clean.index[2]-Hm0_Te_clean.index[1]).seconds<jupyter_output><empty_output><jupyter_text>1. SamplingThe first step is sampling the sea state to get samples $(H_s, T_e)_i$ and associtated weights. For this we will use the `waves.contours.samples_full_seastate` function. We will sample 20 points between each return level, for 10 levels ranging from 0.001—100 years return periods. For more details on the sampling approach see> [1] Coe, <NAME>., <NAME>, <NAME>, and <NAME>. 2018. “Full Long-Term Design Response Analysis of a Wave Energy Converter.” Renewable Energy 116: 356–66.> [2] Eckert-Gallup, <NAME>., <NAME>, <NAME>, and <NAME>. 2016. “Application of Principal Component Analysis (PCA) and Improved Joint Probability Distributions to the Inverse First-Order Reliability Method (I-FORM) for Predicting Extreme Sea States.” Ocean Engineering 112 (January): 307–19.<jupyter_code># return levels
levels = np.array([0.001, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100])

# points per return level interval
npoints = 20

# Create samples
sample_hs, sample_te, sample_weights = contours.samples_full_seastate(
    Hm0, Te, npoints, levels, dt)<jupyter_output><empty_output><jupyter_text>We will now plot the samples alongside the contours. First we will create the different contours using `contours.environmental_contours`. See `environmental_contours_example.ipynb` for more details on using this function. There are 20 samples, randomly distributed, between each set of return levels.<jupyter_code># Create the contours
Te_contours = []
Hm0_contours = []

for period in levels:
    copula = contours.environmental_contours(
        Hm0, Te, dt, period, 'PCA', return_PCA=True)
    Hm0_contours.append(copula['PCA_x1'])
    Te_contours.append(copula['PCA_x2'])

# plot
fig, ax = plt.subplots(figsize=(8, 4))

labels = [f"{period}-year Contour" for period in levels]

ax = graphics.plot_environmental_contour(
    sample_te, sample_hs, Te_contours, Hm0_contours,
    data_label='Samples', contour_label=labels,
    x_label='Energy Period, $Te$ [s]',
    y_label='Sig. wave height, $Hm0$ [m]', ax=ax)<jupyter_output><empty_output><jupyter_text>2. Short-Term Extreme DistributionsMany different methods for short-term extremes were adapted from WDRT, and a summary and examples can be found in `short_term_extremes_example.ipynb`. The response quantity of interest is typically related to the WEC itself, e.g. maximum heave displacement, PTO extension, or load on the mooring lines. This requires running a simulation (e.g. WEC-Sim) for each of the 200 sampled sea states $(H_s, T_e)_i$. For the sake of example we will consider the wave elevation as the quantity of interest (can be thought as a proxy for heave motion in this example). Wave elevation time-series for a specific sea state can be created quickly without running any external software. **NOTE:** The majority of the for loop below is simply creating the synthetic data (wave elevation time series). In a realistic case the variables `time` and `data` describing each time series would be obtained externally, e.g. through simulation software such as WEC-Sim or CFD. For this reason the details of creating the synthetic data are not presented here, instead assume for each sea state there is time-series data available. The last lines of the for-loop create the short-term extreme distribution from the time-series using the `loads.extreme.short_term_extreme` function. The short-term period will be 3-hours and we will use 1-hour "simulations" and the Weibul-tail-fitting method to estimate the 3-hour short-term extreme distributions for each of the 200 samples.For more details on short-term extreme distributions see `short_term_extremes_example.ipynb` and > [3] <NAME>, <NAME>., and <NAME>. 2015. “Comparison of Methods for Estimating Short-Term Extreme Response of Wave Energy Converters.” In OCEANS 2015 - MTS/IEEE Washington, 1–6. IEEE.<jupyter_code># create the short-term extreme distribution for each sample sea state
t_st = 3.0 * 60.0 * 60.0
gamma = 3.3 
t_sim = 1.0 * 60.0 * 60.0

ste_all = []
i = 0
n = len(sample_hs)
for hs, te in zip(sample_hs, sample_te):
    tp = te / (0.8255 + 0.03852*gamma - 0.005537*gamma**2 + 0.0003154*gamma**3)
    i += 1
    print(f"Sea state {i}/{n}. (Hs, Te) = ({hs} m, {te} s). Tp = {tp} s")
    # time & frequency arrays
    f0 = 1.0/t_sim
    T_min = tp/10.0  # s
    f_max = 1.0/T_min
    Nf = int(f_max/f0)
    time = np.linspace(0, t_sim, 2*Nf+1)
    f = np.linspace(f0, f_max, Nf)
    # spectrum
    S = resource.jonswap_spectrum(f, tp, hs, gamma)
    # 1-hour elevation time-series
    data = resource.surface_elevation(S, time).values.squeeze()
    # 3-hour extreme distribution
    ste = extreme.short_term_extreme(time, data, t_st, 'peaks_weibull_tail_fit')
    ste_all.append(ste)<jupyter_output>Sea state 1/200. (Hs, Te) = (2.7476690431712387 m, 9.898396804782426 s). Tp = 10.953763434059923 s
Sea state 2/200. (Hs, Te) = (3.531570115440943 m, 11.072854690290276 s). Tp = 12.253441967345596 s
Sea state 3/200. (Hs, Te) = (3.5631169008261416 m, 10.529043453084103 s). Tp = 11.651649600094585 s
Sea state 4/200. (Hs, Te) = (3.3136880458617415 m, 9.477416257821945 s). Tp = 10.48789795981279 s
Sea state 5/200. (Hs, Te) = (3.172687890319087 m, 9.360633157328573 s). Tp = 10.35866345031299 s
Sea state 6/200. (Hs, Te) = (3.0630002153170452 m, 8.59765613746172 s). Tp = 9.514337854353085 s
Sea state 7/200. (Hs, Te) = (2.7069002802048345 m, 8.619396416083681 s). Tp = 9.538396080519615 s
Sea state 8/200. (Hs, Te) = (2.4473757676516734 m, 8.436380373858702 s). Tp = 9.335866875971888 s
Sea state 9/200. (Hs, Te) = (2.046394958916393 m, 7.376850430979894 s). Tp = 8.163369897472288 s
Sea state 10/200. (Hs, Te) = (2.293289736769183 m, 8.853033813999282 s). Tp = 9.796943887461476 s
Sea state 11/200. ([...]<jupyter_text>3. Long-Term Extreme DistributionFinally we integrate the weighted short-term extreme distributions over the entire sea state space to obtain the extreme distribution, assuming a 3-hour sea state coherence. For this we use the `loads.extreme.full_seastate_long_term_extreme` function.The integral reduces to a sum over the 200 bins of the weighted short-term extreme distributions.<jupyter_code>lte = extreme.full_seastate_long_term_extreme(ste_all, sample_weights)<jupyter_output><empty_output><jupyter_text>Similar to the short-term extreme functions, the output of long-term extreme function is a probability distribution (`scipy.stats.rv_continuous`). This object provides common statistical functions (PDF, CDF, PPF, etc.) and metrics (expected value, median, etc). Here, we will look at the survival function and the 100-year return level. The value of the survival function at a given return level (e.g. 100-years) (`s_t` in the code below) is $1/N$ where $N$ is the number of short-term periods in the return period. In this case $N$ is the number of 3-hour periods in a 100-year period, which gives `s_t`$\approx 3e-6$. The corresponding response, i.e. the 100-year wave elevation in this case, is given as the inverse cumulative function (ppf) of $1-$`s_t`. This gives a 100-year wave of about 10.4 meters.<jupyter_code>t_st_hr = t_st/(60.0*60.0)
t_return_yr = 100.0
s_t = 1.0/(365.25*24*t_return_yr/t_st_hr)
x_t = lte.ppf(1-s_t)

print(f"100-year elevation: {x_t} m")<jupyter_output>100-year elevation: 11.139864671857818 m<jupyter_text>Finally we plot the survival function and show the 100-year return level (dashed grey lines). The 100-year value is about 104 m (where the grey line intersects the x-axis)<jupyter_code>x = np.linspace(0, 20, 1000)

fig, ax = plt.subplots()

# plot survival function
ax.semilogy(x, lte.sf(x))

# format plot
plt.grid(True, which="major", linestyle=":")
ax.tick_params(axis="both", which="major", direction="in")
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both') 
plt.minorticks_off()
ax.set_xticks([0, 5, 10, 15, 20])
ax.set_yticks(1.0*10.0**(-1*np.arange(11)))
ax.set_xlabel("elevation [m]")
ax.set_ylabel("survival function (1-cdf)")
ax.set_xlim([0, x[-1]])
ylim = [1e-10, 1]
ax.set_ylim(ylim)

# 100-year return level
ax.plot([0, x[-1]], [s_t, s_t], '--', color="0.5", linewidth=1)
ax.plot([x_t, x_t], ylim, '--', color="0.5", linewidth=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10142942737310623, 'fraction_numerical': 0.07455276898057006, 'mean_word_length': 4.889112903225806, 'pattern_counts': {'":': 1, '<': 48, '<?xml version=': 0, '>': 50, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '518882', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4523, 'n_tokens_neox': 3882, 'n_words': 1367, 'path': 'examples/extreme_response_full_sea_state_example.ipynb', 'repo_name': 'MHKiT-Code-Hub/mhkit-python', 'size': '11821'}