diff --git a/OSWEC_Optimization_Damping/.gitignore b/OSWEC_Optimization_Damping/.gitignore new file mode 100644 index 00000000..2a56ff95 --- /dev/null +++ b/OSWEC_Optimization_Damping/.gitignore @@ -0,0 +1,51 @@ +# MATLAB autosave/backup +*.asv +*.m~ +*.bak + +# MATLAB / Simulink generated files +*.slxc +slprj/ +codegen/ +work/ +*.log +*.fig + +# Temporary batch files +currentRun.mat +latest_run.mat +latest_damping_bounds.mat +oswec_batch.mat +oswec_results.mat +oswec_optimization_summary.mat + +# Generated run outputs +runs/* +!runs/.gitkeep + +# Generated damping-bound outputs +damping_bounds/* +!damping_bounds/.gitkeep + +# Generated result files +oswec_*_optimization_summary.mat +oswec_*_optimization_summary.csv +oswec_*_full_run_results.csv +oswec_*_best_weighted_condition_results.csv +oswec_*_best_damping_by_condition.csv +oswec_*_weighted_vs_unweighted_power.png +oswec_*_failed_runs_by_damping.png +oswec_*_power_by_condition.png +oswec_*_weighted_contribution_best_damping.png + +# Python cache +__pycache__/ +*.pyc +.ipynb_checkpoints/ + +# Generated wave-condition files +wave_conditions/*.mat + +# OS files +.DS_Store +Thumbs.db diff --git a/OSWEC_Optimization_Damping/OSWEC.slx b/OSWEC_Optimization_Damping/OSWEC.slx new file mode 100644 index 00000000..42d0e230 Binary files /dev/null and b/OSWEC_Optimization_Damping/OSWEC.slx differ diff --git a/OSWEC_Optimization_Damping/README.md b/OSWEC_Optimization_Damping/README.md new file mode 100644 index 00000000..48687aeb --- /dev/null +++ b/OSWEC_Optimization_Damping/README.md @@ -0,0 +1,21 @@ +**Author:** Rebekah A. Saucier + +**Geometry:** OSWEC + +**Original Version:** WEC-Sim v6.1 + +**Dependencies:** +* WEC-Sim +* MATLAB/Simulink +* Precomputed MHKiT wave-condition `.mat` files in `wave_conditions/` + +**Description** + +The **OSWEC_Optimization_Damping** example uses WEC-Sim multiple condition runs, MCR, to sweep PTO damping values for the OSWEC model. Representative wave conditions are precomputed from MHKiT/NDBC wave resource data and stored in `wave_conditions/`. + +The input file expands the selected wave-condition file across user-defined PTO damping values. Running `wecSimMCR` simulates all wave-condition and damping combinations. The `userDefinedFunctions.m` file suppresses per-case plots during MCR runs and summarizes weighted mean power, unweighted mean power, AEP, and the best damping value. + +To change the study, edit the selected wave-condition file and `dampingValues` in `wecSimInputFile.m`, then run: + +```matlab +wecSimMCR diff --git a/OSWEC_Optimization_Damping/Waveconditions_WEC-Sim.ipynb b/OSWEC_Optimization_Damping/Waveconditions_WEC-Sim.ipynb new file mode 100644 index 00000000..008a55ed --- /dev/null +++ b/OSWEC_Optimization_Damping/Waveconditions_WEC-Sim.ipynb @@ -0,0 +1,655 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "51d7c10e-c5b9-430a-88c3-c332140b31c3", + "metadata": {}, + "source": [ + "## Functionality of this script\n", + " 1. Downloads NDBC spectral wave density data\n", + " 2. Computes Hm0, Te, Tp, J, and Tz using MHKiT\n", + " 3. Clusters the sea states into X representative wave conditions\n", + " 4. Saves precomputed `.mat` and `.csv` files for WEC-Sim/MATLAB\n", + "## For WEC-Sim damping optimization, the most important saved variables are from this framework are:\n", + " - condition_id\n", + " - H = Hm0\n", + " - T = Tp\n", + " - Hm0\n", + " - Te\n", + " - Tp\n", + " - weights\n", + " - probability\n", + "For full documentation please see the MHKIT PacWave Assesment here:\n", + "\n", + "https://mhkit-software.github.io/MHKiT/PacWave_resource_characterization_example.html" + ] + }, + { + "cell_type": "markdown", + "id": "97b08b10-2647-41ef-b1f7-9fed864a1fba", + "metadata": {}, + "source": [ + "Imported package dependencies needed for this script." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c4bde20f-31f6-4f52-a8e5-03129a0abd80", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import os\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.mixture import GaussianMixture\n", + "from scipy.io import savemat\n", + "\n", + "from mhkit.wave import resource\n", + "from mhkit.wave.io import ndbc" + ] + }, + { + "cell_type": "markdown", + "id": "7290195f-6dc9-4fa8-828c-561ef100b86d", + "metadata": {}, + "source": [ + "User inputs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08bfd705-1d75-4c55-889d-774680225043", + "metadata": {}, + "outputs": [], + "source": [ + "# Change this to your project location\n", + "PROJECT_ROOT = Path.home() / \"WEC-Sim_Applications\" / \"OSWEC_Optimization_Damping\"\n", + "\n", + "# NDBC buoy number\n", + "buoy_number = \"46050\"\n", + "\n", + "# Years to analyze\n", + "years = [\"2020\", \"2021\", \"2022\", \"2023\", \"2024\", \"2025\"]\n", + "\n", + "# Number of representative wave conditions to create.\n", + "# These are your X values.\n", + "clusters = [4, 8, 16, 32, 64]\n", + "\n", + "# Water depth in meters\n", + "water_depth = 160.0\n", + "\n", + "# Choose which cluster result to export.\n", + "# Set to a number like 32 to export only one file.\n", + "# Set to None to export all cluster cases listed above.\n", + "cluster_to_export = None\n", + "\n", + "# Make diagnostic cluster plots?\n", + "# Keep False for clean production runs.\n", + "make_plots = False\n", + "\n", + "# Random seed for reproducible clustering\n", + "random_state = 1\n" + ] + }, + { + "cell_type": "markdown", + "id": "92b3c604-410b-4a01-bccb-d7c8888f7b3e", + "metadata": {}, + "source": [ + "Project folder set up" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38111987-c50d-444b-a51e-f5fdbff11976", + "metadata": {}, + "outputs": [], + "source": [ + "os.chdir(PROJECT_ROOT)\n", + "\n", + "output_folder = PROJECT_ROOT / \"wave_conditions\"\n", + "output_folder.mkdir(parents=True, exist_ok=True)\n", + "\n", + "print(\"--------------------------------------------\")\n", + "print(\"Current working directory:\")\n", + "print(Path.cwd())\n", + "print(\"--------------------------------------------\")\n", + "print(f\"Project root: {PROJECT_ROOT}\")\n", + "print(f\"Output folder: {output_folder}\")\n", + "print(\"--------------------------------------------\")" + ] + }, + { + "cell_type": "markdown", + "id": "dd5e1448-7ced-4ef9-9baa-32365fbcf854", + "metadata": {}, + "source": [ + "Download the spectral wave density data from the NDBC" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "588dd552-d407-4585-acbe-9870dea70a99", + "metadata": {}, + "outputs": [], + "source": [ + "parameter = \"swden\"\n", + "\n", + "print(f\"Buoy number: {buoy_number}\")\n", + "print(f\"Years: {years}\")\n", + "print(f\"Clusters: {clusters}\")\n", + "print(f\"Water depth: {water_depth} m\")\n", + "print(\"--------------------------------------------\")\n", + "\n", + "print(\"Checking available NDBC data...\")\n", + "\n", + "ndbc_available_data = ndbc.available_data(parameter, buoy_number)\n", + "\n", + "# Create a clean year string for filenames\n", + "if years is None:\n", + " years_string = \"all_years\"\n", + "else:\n", + " years = [str(year) for year in years]\n", + " years_string = \"_\".join(years)\n", + "\n", + " # Filter available files by selected years\n", + " year_pattern = \"|\".join(years)\n", + "\n", + " ndbc_available_data = ndbc_available_data[\n", + " ndbc_available_data[\"filename\"].astype(str).str.contains(year_pattern)\n", + " ]\n", + "\n", + " if ndbc_available_data.empty:\n", + " raise ValueError(\n", + " f\"No data files found for buoy {buoy_number} and years {years}.\"\n", + " )\n", + "\n", + "filenames = ndbc_available_data[\"filename\"]\n", + "\n", + "print(\"Downloading NDBC data...\")\n", + "\n", + "ndbc_requested_data = ndbc.request_data(parameter, filenames)" + ] + }, + { + "cell_type": "markdown", + "id": "bd5ebb06-11f2-4125-b17f-daeb738d7dc3", + "metadata": {}, + "source": [ + "Clean NDBC dara and create datetime index" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d94b1d9f-e555-47ca-89a6-4b2b94377269", + "metadata": {}, + "outputs": [], + "source": [ + "ndbc_data = {}\n", + "\n", + "for year in ndbc_requested_data:\n", + " print(f\"Cleaning data for year {year}...\")\n", + "\n", + " data_raw = ndbc_requested_data[year].copy()\n", + "\n", + " # Create datetime index from NDBC date/time columns\n", + " data_raw[\"date\"] = pd.to_datetime(\n", + " {\n", + " \"year\": data_raw[\"#YY\"],\n", + " \"month\": data_raw[\"MM\"],\n", + " \"day\": data_raw[\"DD\"],\n", + " \"hour\": data_raw[\"hh\"],\n", + " \"minute\": data_raw[\"mm\"],\n", + " },\n", + " errors=\"coerce\",\n", + " )\n", + "\n", + " data_raw = data_raw.set_index(\"date\")\n", + "\n", + " # Drop original date/time columns\n", + " data_raw = data_raw.drop(columns=[\"#YY\", \"MM\", \"DD\", \"hh\", \"mm\"])\n", + "\n", + " # Convert frequency column names to floats\n", + " new_columns = []\n", + "\n", + " for col in data_raw.columns:\n", + " col_string = str(col)\n", + "\n", + " if col_string.startswith(\".\"):\n", + " col_string = \"0\" + col_string\n", + "\n", + " new_columns.append(float(col_string))\n", + "\n", + " data_raw.columns = new_columns\n", + "\n", + " # Convert data values to numeric\n", + " data_raw = data_raw.apply(pd.to_numeric, errors=\"coerce\")\n", + "\n", + " # Replace NDBC missing/bad values\n", + " data_raw = data_raw.replace([999.0, 99.0], np.nan)\n", + "\n", + " # Drop rows with missing data\n", + " data_raw = data_raw.dropna()\n", + "\n", + " # Sort by time\n", + " data_raw = data_raw.sort_index()\n", + "\n", + " ndbc_data[str(year)] = data_raw\n", + "\n", + "print(\"Cleaned data years:\")\n", + "print(list(ndbc_data.keys()))" + ] + }, + { + "cell_type": "markdown", + "id": "8c03146a-b1d0-499e-8c0c-35afcaa76903", + "metadata": {}, + "source": [ + "Calculate Quantities of Interest (QOI's)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "292cd10d-756e-494a-95e5-b7cfe0f41e37", + "metadata": {}, + "outputs": [], + "source": [ + "Hm0_list = []\n", + "Te_list = []\n", + "J_list = []\n", + "Tp_list = []\n", + "Tz_list = []\n", + "\n", + "for year in ndbc_data:\n", + " print(f\"Calculating QoIs for year {year}...\")\n", + "\n", + " data_raw = ndbc_data[year]\n", + "\n", + " year_data = data_raw[data_raw != 999.0].dropna()\n", + "\n", + " # MHKiT expects frequency as rows and timestamps as columns\n", + " spectrum = year_data.T\n", + "\n", + " Hm0_list.append(resource.significant_wave_height(spectrum))\n", + " Te_list.append(resource.energy_period(spectrum))\n", + " J_list.append(resource.energy_flux(spectrum, h=water_depth))\n", + "\n", + " # Peak period calculation from spectral peak frequency\n", + " fp = spectrum.idxmax(axis=0).astype(float)\n", + " Tp = 1.0 / fp\n", + " Tp = pd.DataFrame(Tp, index=spectrum.columns, columns=[\"Tp\"])\n", + " Tp_list.append(Tp)\n", + "\n", + " Tz_list.append(resource.average_zero_crossing_period(spectrum))\n", + "\n", + "\n", + "Te = pd.concat(Te_list, axis=0)\n", + "Tp = pd.concat(Tp_list, axis=0)\n", + "Hm0 = pd.concat(Hm0_list, axis=0)\n", + "J = pd.concat(J_list, axis=0)\n", + "Tz = pd.concat(Tz_list, axis=0)\n", + "\n", + "# Combine into one DataFrame\n", + "data = pd.concat([Hm0, Te, Tp, J, Tz], axis=1)\n", + "\n", + "# Make sure columns are named correctly\n", + "data.columns = [\"Hm0\", \"Te\", \"Tp\", \"J\", \"Tz\"]\n", + "\n", + "# Calculate mean wave steepness\n", + "data[\"Sm\"] = data.Hm0 / (9.81 / (2.0 * np.pi) * data.Tz**2)\n", + "\n", + "# Drop NaNs and sort\n", + "data.dropna(inplace=True)\n", + "data.sort_index(inplace=True)\n", + "\n", + "print(\"--------------------------------------------\")\n", + "print(\"QoI data preview:\")\n", + "print(data.head())\n", + "print(\"--------------------------------------------\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "82aa954c-af26-43f1-b481-34afef1ede84", + "metadata": {}, + "source": [ + "Clean up " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bdc022c8-033a-471b-8668-6b7d9836dce2", + "metadata": {}, + "outputs": [], + "source": [ + "data_clean = data.copy()\n", + "\n", + "# Remove unrealistic wave heights\n", + "data_clean = data_clean[data_clean.Hm0 < 20]\n", + "\n", + "# Keep your original cleaning approach\n", + "sigma = data_clean.J.std()\n", + "data_clean = data_clean[data_clean.J > (data_clean.J.mean() - 0.9 * sigma)]\n", + "\n", + "print(f\"Number of sea states after cleaning: {len(data_clean)}\")\n", + "print(\"--------------------------------------------\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "3b071219-49d6-4f5c-b278-3df14391b536", + "metadata": {}, + "source": [ + "Pacwave style Clustering Using Gaussian mixture model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa394cf0-63f7-42cb-874c-761dbc36ae89", + "metadata": {}, + "outputs": [], + "source": [ + "# Cluster on energy period and significant wave height\n", + "X = np.vstack((data_clean.Te.values, data_clean.Hm0.values)).T\n", + "\n", + "results = {}\n", + "\n", + "if make_plots:\n", + " fig, axs = plt.subplots(\n", + " len(clusters),\n", + " 1,\n", + " figsize=(8, 5 * len(clusters)),\n", + " sharex=True,\n", + " )\n", + "\n", + " if len(clusters) == 1:\n", + " axs = [axs]\n", + "\n", + "\n", + "for plot_index, cluster in enumerate(clusters):\n", + " print(f\"Creating {cluster} representative wave conditions...\")\n", + "\n", + " gmm = GaussianMixture(\n", + " n_components=cluster,\n", + " random_state=random_state,\n", + " ).fit(X)\n", + "\n", + " labels = gmm.predict(X)\n", + "\n", + " # Save cluster centers and weights\n", + " result = pd.DataFrame(gmm.means_, columns=[\"Te\", \"Hm0\"])\n", + " result[\"weights\"] = gmm.weights_\n", + "\n", + " # Normalize weights just to be safe\n", + " result[\"weights\"] = result[\"weights\"] / result[\"weights\"].sum()\n", + "\n", + " # Use the same relationship as the PacWave example\n", + " result[\"Tp\"] = result.Te / 0.858\n", + "\n", + " # Sort from smaller waves to larger waves for easier interpretation\n", + " result = result.sort_values([\"Hm0\", \"Te\"]).reset_index(drop=True)\n", + "\n", + " # Add condition IDs after sorting\n", + " result.insert(0, \"condition_id\", np.arange(1, len(result) + 1))\n", + "\n", + " # Add WEC-Sim-friendly aliases\n", + " result[\"H\"] = result[\"Hm0\"]\n", + " result[\"T\"] = result[\"Tp\"]\n", + "\n", + " # Use probability as clearer name for optimization weighting\n", + " result[\"probability\"] = result[\"weights\"]\n", + "\n", + " results[cluster] = result\n", + "\n", + " if make_plots:\n", + " axs[plot_index].scatter(\n", + " data_clean.Te.values,\n", + " data_clean.Hm0.values,\n", + " c=labels,\n", + " s=8,\n", + " )\n", + "\n", + " axs[plot_index].plot(\n", + " result.Te,\n", + " result.Hm0,\n", + " \"m+\",\n", + " markersize=10,\n", + " )\n", + "\n", + " axs[plot_index].set_title(f\"{cluster} Clusters\")\n", + " axs[plot_index].set_ylabel(\"Hm0 [m]\")\n", + "\n", + "if make_plots:\n", + " axs[-1].set_xlabel(\"Energy Period, Te [s]\")\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d03b41ce-a0fa-4f14-ba9d-5f1d35d656f0", + "metadata": {}, + "source": [ + "View selected wave conditions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95d6601a-44e6-4a40-bc89-6c95e65688ac", + "metadata": {}, + "outputs": [], + "source": [ + "if cluster_to_export is not None:\n", + " print(\"--------------------------------------------\")\n", + " print(f\"Representative wave conditions for {cluster_to_export} clusters:\")\n", + " print(\"--------------------------------------------\")\n", + " print(results[cluster_to_export])" + ] + }, + { + "cell_type": "markdown", + "id": "47b159fc-de2a-4bb9-84b2-6974a49c7a9f", + "metadata": {}, + "source": [ + "Save the representative wave conditions as `.mat` and `.csv` file types " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c0d99572-8a63-4d8c-bbd9-32910bca1b12", + "metadata": {}, + "outputs": [], + "source": [ + "def make_matlab_mcr_struct(result):\n", + " \"\"\"\n", + " Create a MATLAB-compatible struct array named mcr.\n", + "\n", + " Each mcr(i) contains one wave condition:\n", + " mcr(i).caseNum\n", + " mcr(i).condition_id\n", + " mcr(i).H\n", + " mcr(i).T\n", + " mcr(i).Hm0\n", + " mcr(i).Te\n", + " mcr(i).Tp\n", + " mcr(i).weights\n", + " mcr(i).probability\n", + " mcr(i).caseName\n", + " \"\"\"\n", + "\n", + " n_cases = len(result)\n", + "\n", + " mcr_dtype = np.dtype(\n", + " [\n", + " (\"caseNum\", \"O\"),\n", + " (\"condition_id\", \"O\"),\n", + " (\"H\", \"O\"),\n", + " (\"T\", \"O\"),\n", + " (\"Hm0\", \"O\"),\n", + " (\"Te\", \"O\"),\n", + " (\"Tp\", \"O\"),\n", + " (\"weights\", \"O\"),\n", + " (\"probability\", \"O\"),\n", + " (\"caseName\", \"O\"),\n", + " ]\n", + " )\n", + "\n", + " # 1-by-n MATLAB struct array\n", + " mcr = np.empty((1, n_cases), dtype=mcr_dtype)\n", + "\n", + " for i in range(n_cases):\n", + " mcr[\"caseNum\"][0, i] = i + 1\n", + " mcr[\"condition_id\"][0, i] = float(result[\"condition_id\"].iloc[i])\n", + "\n", + " # WEC-Sim wave inputs\n", + " mcr[\"H\"][0, i] = float(result[\"H\"].iloc[i])\n", + " mcr[\"T\"][0, i] = float(result[\"T\"].iloc[i])\n", + "\n", + " # Resource variables\n", + " mcr[\"Hm0\"][0, i] = float(result[\"Hm0\"].iloc[i])\n", + " mcr[\"Te\"][0, i] = float(result[\"Te\"].iloc[i])\n", + " mcr[\"Tp\"][0, i] = float(result[\"Tp\"].iloc[i])\n", + "\n", + " # Occurrence weighting\n", + " mcr[\"weights\"][0, i] = float(result[\"weights\"].iloc[i])\n", + " mcr[\"probability\"][0, i] = float(result[\"probability\"].iloc[i])\n", + "\n", + " # Human-readable case name\n", + " mcr[\"caseName\"][0, i] = f\"condition_{int(result['condition_id'].iloc[i]):03d}\"\n", + "\n", + " return mcr\n", + "\n", + "\n", + "# Decide which cluster counts to export\n", + "if cluster_to_export is None:\n", + " clusters_to_save = clusters\n", + "else:\n", + " clusters_to_save = [cluster_to_export]\n", + "\n", + "for cluster in clusters_to_save:\n", + " result = results[cluster]\n", + "\n", + " # MATLAB/WEC-Sim-friendly array output\n", + " condition_id = result[\"condition_id\"].to_numpy(dtype=float)\n", + "\n", + " Hm0_array = result[\"Hm0\"].to_numpy(dtype=float)\n", + " Te_array = result[\"Te\"].to_numpy(dtype=float)\n", + " Tp_array = result[\"Tp\"].to_numpy(dtype=float)\n", + "\n", + " weights_array = result[\"weights\"].to_numpy(dtype=float)\n", + " probability_array = result[\"probability\"].to_numpy(dtype=float)\n", + "\n", + " # WEC-Sim aliases\n", + " H_array = result[\"H\"].to_numpy(dtype=float)\n", + " T_array = result[\"T\"].to_numpy(dtype=float)\n", + "\n", + " # Plain array-style MAT file\n", + " wave_conditions = {\n", + " \"condition_id\": condition_id,\n", + " \"Hm0\": Hm0_array,\n", + " \"Te\": Te_array,\n", + " \"Tp\": Tp_array,\n", + " \"H\": H_array,\n", + " \"T\": T_array,\n", + " \"weights\": weights_array,\n", + " \"probability\": probability_array,\n", + " \"n_conditions\": np.array([cluster], dtype=float),\n", + " \"cluster_to_export\": np.array([cluster], dtype=float),\n", + " \"buoy_number\": buoy_number,\n", + " \"years_string\": years_string,\n", + " \"water_depth\": np.array([water_depth], dtype=float),\n", + " }\n", + "\n", + " # MCR-style MAT file\n", + " mcr = make_matlab_mcr_struct(result)\n", + "\n", + " output_base_name = (\n", + " f\"wave_conditions_buoy_{buoy_number}_{years_string}_{cluster}_clusters\"\n", + " )\n", + "\n", + " output_mat_file = output_folder / f\"{output_base_name}.mat\"\n", + " output_mcr_file = output_folder / f\"{output_base_name}_mcr.mat\"\n", + " output_csv_file = output_folder / f\"{output_base_name}.csv\"\n", + "\n", + " # Save normal MAT file\n", + " savemat(str(output_mat_file), wave_conditions)\n", + "\n", + " # Save WEC-Sim MCR MAT file\n", + " savemat(\n", + " str(output_mcr_file),\n", + " {\n", + " \"mcr\": mcr,\n", + " \"condition_id\": condition_id,\n", + " \"H\": H_array,\n", + " \"T\": T_array,\n", + " \"Hm0\": Hm0_array,\n", + " \"Te\": Te_array,\n", + " \"Tp\": Tp_array,\n", + " \"weights\": weights_array,\n", + " \"probability\": probability_array,\n", + " \"n_conditions\": np.array([cluster], dtype=float),\n", + " \"cluster_to_export\": np.array([cluster], dtype=float),\n", + " \"buoy_number\": buoy_number,\n", + " \"years_string\": years_string,\n", + " \"water_depth\": np.array([water_depth], dtype=float),\n", + " },\n", + " )\n", + "\n", + " # Save CSV for easy inspection\n", + " result.to_csv(output_csv_file, index=False)\n", + "\n", + " print(\"--------------------------------------------\")\n", + " print(f\"Saved {cluster}-condition WEC-Sim wave files:\")\n", + " print(f\"Array MAT: {output_mat_file}\")\n", + " print(f\"MCR MAT: {output_mcr_file}\")\n", + " print(f\"CSV: {output_csv_file}\")\n", + " print(\"--------------------------------------------\")\n", + "\n", + "print(\"Done.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58310ff5-65ed-4049-bf42-82c3e1ebe989", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/OSWEC_Optimization_Damping/elevationData.mat b/OSWEC_Optimization_Damping/elevationData.mat new file mode 100644 index 00000000..ffea9b7c Binary files /dev/null and b/OSWEC_Optimization_Damping/elevationData.mat differ diff --git a/OSWEC_Optimization_Damping/generate_wave_conditions.ipynb b/OSWEC_Optimization_Damping/generate_wave_conditions.ipynb new file mode 100644 index 00000000..e614001f --- /dev/null +++ b/OSWEC_Optimization_Damping/generate_wave_conditions.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "5968d78a-a7f8-4126-af24-038aa471b625", + "metadata": {}, + "outputs": [], + "source": [ + "# Simplified PacWave-style representative wave conditions\n", + "# for WEC-Sim optimization\n", + "# saves .mat files to /wave_conditions\n", + "\n", + "\n", + "from pathlib import Path\n", + "import os\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.mixture import GaussianMixture\n", + "from scipy.io import savemat\n", + "\n", + "from mhkit.wave import resource\n", + "from mhkit.wave.io import ndbc\n", + "# Project folder setup\n", + "#Change the PROJECT Root for your folder location\n", + "PROJECT_ROOT = Path.home() / \"WEC-Sim\" / \"OSWEC_Optimization_Damping\" #CHANGE THIS TO YOUR LOCATION\n", + "# Change Python/Jupyter working directory to project root\n", + "os.chdir(PROJECT_ROOT)\n", + "\n", + "# Save wave-condition files inside the project folder\n", + "output_folder = PROJECT_ROOT / \"wave_conditions\"\n", + "output_folder.mkdir(parents=True, exist_ok=True)\n", + "\n", + "print(\"--------------------------------------------\")\n", + "print(\"Current working directory:\")\n", + "print(Path.cwd())\n", + "print(\"--------------------------------------------\")\n", + "print(f\"Project root: {PROJECT_ROOT}\")\n", + "print(f\"Output folder: {output_folder}\")\n", + "print(\"--------------------------------------------\")\n", + "\n", + "# INPUTS\n", + "# NDBC buoy number\n", + "buoy_number = \"46050\"\n", + "\n", + "# Years to analyze\n", + "# Example: [\"2020\"] or [\"2019\", \"2020\", \"2021\"]\n", + "years = [\"2020\", \"2021\", \"2022\", \"2023\", \"2024\", \"2025\"]\n", + "\n", + "# Number of representative wave conditions to create\n", + "\n", + "clusters = [4, 8, 16, 32, 64]\n", + "\n", + "# Water depth in meters\n", + "# Use the buoy/site water depth you want for energy flux calculation\n", + "#(see MHKIT example for more details on this)\n", + "water_depth = 160.0\n", + "\n", + "# Choose which cluster result to export to .mat for WEC-Sim\n", + "# Set to a number like 32 to export only one file\n", + "# Set to None to export all clusters\n", + "cluster_to_export = None\n", + "# Download spectral wave density data from NDBC\n", + "\n", + "parameter = \"swden\"\n", + "\n", + "print(f\"Buoy number: {buoy_number}\")\n", + "print(f\"Years: {years}\")\n", + "print(f\"Clusters: {clusters}\")\n", + "print(f\"Water depth: {water_depth} m\")\n", + "print(\"--------------------------------------------\")\n", + "\n", + "print(\"Checking available NDBC data...\")\n", + "\n", + "ndbc_available_data = ndbc.available_data(parameter, buoy_number)\n", + "\n", + "# Create a clean year string for filenames\n", + "if years is None:\n", + " years_string = \"all_years\"\n", + "else:\n", + " years = [str(year) for year in years]\n", + " years_string = \"_\".join(years)\n", + "\n", + " # Filter available files by selected years\n", + " year_pattern = \"|\".join(years)\n", + "\n", + " ndbc_available_data = ndbc_available_data[\n", + " ndbc_available_data[\"filename\"].astype(str).str.contains(year_pattern)\n", + " ]\n", + "\n", + " if ndbc_available_data.empty:\n", + " raise ValueError(f\"No data files found for buoy {buoy_number} and years {years}.\")\n", + "\n", + "filenames = ndbc_available_data[\"filename\"]\n", + "\n", + "print(\"Downloading NDBC data...\")\n", + "\n", + "ndbc_requested_data = ndbc.request_data(parameter, filenames)\n", + "\n", + "# Clean NDBC data and create DateTime index\n", + "\n", + "\n", + "ndbc_data = {}\n", + "\n", + "for year in ndbc_requested_data:\n", + " print(f\"Cleaning data for year {year}...\")\n", + "\n", + " data_raw = ndbc_requested_data[year].copy()\n", + "\n", + " # Create datetime index from NDBC date/time columns\n", + " data_raw[\"date\"] = pd.to_datetime(\n", + " {\n", + " \"year\": data_raw[\"#YY\"],\n", + " \"month\": data_raw[\"MM\"],\n", + " \"day\": data_raw[\"DD\"],\n", + " \"hour\": data_raw[\"hh\"],\n", + " \"minute\": data_raw[\"mm\"],\n", + " },\n", + " errors=\"coerce\",\n", + " )\n", + "\n", + " data_raw = data_raw.set_index(\"date\")\n", + "\n", + " # Drop original date/time columns\n", + " data_raw = data_raw.drop(columns=[\"#YY\", \"MM\", \"DD\", \"hh\", \"mm\"])\n", + "\n", + " # Convert frequency column names to floats\n", + " new_columns = []\n", + "\n", + " for col in data_raw.columns:\n", + " col_string = str(col)\n", + "\n", + " if col_string.startswith(\".\"):\n", + " col_string = \"0\" + col_string\n", + "\n", + " new_columns.append(float(col_string))\n", + "\n", + " data_raw.columns = new_columns\n", + "\n", + " # Convert data values to numeric\n", + " data_raw = data_raw.apply(pd.to_numeric, errors=\"coerce\")\n", + "\n", + " # Replace NDBC missing/bad values\n", + " data_raw = data_raw.replace([999.0, 99.0], np.nan)\n", + "\n", + " # Drop rows with missing data\n", + " data_raw = data_raw.dropna()\n", + "\n", + " # Sort by time\n", + " data_raw = data_raw.sort_index()\n", + "\n", + " ndbc_data[str(year)] = data_raw\n", + "\n", + "print(\"Cleaned data years:\")\n", + "print(list(ndbc_data.keys()))\n", + "\n", + "# Calculate QoIs from this data\n", + "\n", + "\n", + "Hm0_list = []\n", + "Te_list = []\n", + "J_list = []\n", + "Tp_list = []\n", + "Tz_list = []\n", + "\n", + "for year in ndbc_data:\n", + " print(f\"Calculating QoIs for year {year}...\")\n", + "\n", + " data_raw = ndbc_data[year]\n", + "\n", + " year_data = data_raw[data_raw != 999.0].dropna()\n", + "\n", + " # MHKiT expects frequency as rows and timestamps as columns\n", + " spectrum = year_data.T\n", + "\n", + " Hm0_list.append(resource.significant_wave_height(spectrum))\n", + " Te_list.append(resource.energy_period(spectrum))\n", + " J_list.append(resource.energy_flux(spectrum, h=water_depth))\n", + "\n", + " # Peak period calculation from PacWave example\n", + " fp = spectrum.idxmax(axis=0).astype(float)\n", + " Tp = 1 / fp\n", + " Tp = pd.DataFrame(Tp, index=spectrum.columns, columns=[\"Tp\"])\n", + " Tp_list.append(Tp)\n", + "\n", + " Tz_list.append(resource.average_zero_crossing_period(spectrum))\n", + "\n", + "\n", + "Te = pd.concat(Te_list, axis=0)\n", + "Tp = pd.concat(Tp_list, axis=0)\n", + "Hm0 = pd.concat(Hm0_list, axis=0)\n", + "J = pd.concat(J_list, axis=0)\n", + "Tz = pd.concat(Tz_list, axis=0)\n", + "\n", + "# Name each Series/DataFrame\n", + "Te.name = \"Te\"\n", + "Tp.name = \"Tp\"\n", + "Hm0.name = \"Hm0\"\n", + "J.name = \"J\"\n", + "Tz.name = \"Tz\"\n", + "\n", + "# Combine into one DataFrame\n", + "data = pd.concat([Hm0, Te, Tp, J, Tz], axis=1)\n", + "\n", + "# Make sure columns are named correctly\n", + "data.columns = [\"Hm0\", \"Te\", \"Tp\", \"J\", \"Tz\"]\n", + "\n", + "# Calculate wave steepness\n", + "data[\"Sm\"] = data.Hm0 / (9.81 / (2 * np.pi) * data.Tz**2)\n", + "\n", + "# Drop NaNs and sort\n", + "data.dropna(inplace=True)\n", + "data.sort_index(inplace=True)\n", + "\n", + "print(\"--------------------------------------------\")\n", + "print(\"QoI data preview:\")\n", + "print(data.head())\n", + "print(\"--------------------------------------------\")\n", + "\n", + "\n", + "data_clean = data[data.Hm0 < 20]\n", + "\n", + "sigma = data_clean.J.std()\n", + "data_clean = data_clean[data_clean.J > (data_clean.J.mean() - 0.9 * sigma)]\n", + "\n", + "print(f\"Number of sea states after cleaning: {len(data_clean)}\")\n", + "\n", + "\n", + "\n", + "# PacWave-style clustering using Gaussian Mixture model\n", + "\n", + "\n", + "X = np.vstack((data_clean.Te.values, data_clean.Hm0.values)).T\n", + "\n", + "fig, axs = plt.subplots(len(clusters), 1, figsize=(8, 5 * len(clusters)), sharex=True)\n", + "\n", + "if len(clusters) == 1:\n", + " axs = [axs]\n", + "\n", + "results = {}\n", + "\n", + "for cluster in clusters:\n", + " gmm = GaussianMixture(n_components=cluster).fit(X)\n", + "\n", + " # Save centers and weights\n", + " result = pd.DataFrame(gmm.means_, columns=[\"Te\", \"Hm0\"])\n", + " result[\"weights\"] = gmm.weights_\n", + "\n", + " # Same relationship used in the PacWave example\n", + " result[\"Tp\"] = result.Te / 0.858\n", + "\n", + " results[cluster] = result\n", + "\n", + " labels = gmm.predict(X)\n", + "\n", + " i = clusters.index(cluster)\n", + "\n", + " axs[i].scatter(\n", + " data_clean.Te.values,\n", + " data_clean.Hm0.values,\n", + " c=labels,\n", + " s=40,\n", + " )\n", + "\n", + " axs[i].plot(\n", + " result.Te,\n", + " result.Hm0,\n", + " \"m+\",\n", + " )\n", + "\n", + " axs[i].title.set_text(f\"{cluster} Clusters\")\n", + " plt.setp(axs[i], ylabel=\"Sig. wave height, $Hm0$ [m]\")\n", + "\n", + "plt.setp(axs[len(clusters) - 1], xlabel=\"Energy Period, $T_e$ [s]\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "\n", + "\n", + "# View the representative wave conditions\n", + "\n", + "\n", + "if cluster_to_export is not None:\n", + " print(\"--------------------------------------------\")\n", + " print(f\"Representative wave conditions for {cluster_to_export} clusters:\")\n", + " print(\"--------------------------------------------\")\n", + " print(results[cluster_to_export])\n", + "\n", + "\n", + "\n", + "# Save representative wave conditions to .mat for WEC-Sim\n", + "\n", + "\n", + "# Decide which cluster counts to export\n", + "if cluster_to_export is None:\n", + " clusters_to_save = clusters\n", + "else:\n", + " clusters_to_save = [cluster_to_export]\n", + "\n", + "for cluster in clusters_to_save:\n", + " result = results[cluster]\n", + "\n", + " wave_conditions = {\n", + " \"Hm0\": result[\"Hm0\"].to_numpy(),\n", + " \"Te\": result[\"Te\"].to_numpy(),\n", + " \"Tp\": result[\"Tp\"].to_numpy(),\n", + " \"weights\": result[\"weights\"].to_numpy(),\n", + " \"cluster_to_export\": np.array([cluster]),\n", + " \"buoy_number\": buoy_number,\n", + " \"years_string\": years_string,\n", + " }\n", + "\n", + " output_mat_file = output_folder / (\n", + " f\"wave_conditions_buoy_{buoy_number}_{years_string}_{cluster}_clusters.mat\"\n", + " )\n", + "\n", + " savemat(str(output_mat_file), wave_conditions)\n", + "\n", + " print(\"--------------------------------------------\")\n", + " print(\"Saved WEC-Sim wave condition file:\")\n", + " print(output_mat_file)\n", + " print(\"--------------------------------------------\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "640cbe04-dea5-4895-b07c-e3c67007c627", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/OSWEC_Optimization_Damping/oswecBuildDampingMCR.m b/OSWEC_Optimization_Damping/oswecBuildDampingMCR.m new file mode 100644 index 00000000..f11e7b30 --- /dev/null +++ b/OSWEC_Optimization_Damping/oswecBuildDampingMCR.m @@ -0,0 +1,165 @@ +function runtimeMcrFile = oswecBuildDampingMCR(baseWaveMcrFile, dampingValues, runtimeMcrFile) +%OSWECBUILDDAMPINGMCR Build WEC-Sim-compatible MCR cases for OSWEC damping sweep. +% +% This function reads a Python/MHKiT-generated wave-condition MAT file and +% expands it into a WEC-Sim MCR file containing: +% +% mcr.header +% mcr.cases +% +% Required wave information: +% Hm0 +% Tp +% +% Optional: +% Te +% probability or weights +% condition_id +% +% The generated mcr.cases columns are: +% 1 caseNum +% 2 condition_id +% 3 wave_id +% 4 damping_id +% 5 Hm0 +% 6 Tp +% 7 Te +% 8 probability +% 9 damping +% +% Example: +% oswecBuildDampingMCR( ... +% fullfile('wave_conditions','wave_conditions_buoy_46050_..._4_clusters.mat'), ... +% [1e3 3e3 1e4 3e4 1e5].', ... +% fullfile('wave_conditions','oswec_runtime_damping_grid_mcr.mat')); + + arguments + baseWaveMcrFile char + dampingValues (:,1) double {mustBePositive} + runtimeMcrFile char + end + + if ~exist(baseWaveMcrFile, 'file') + error('Base wave-condition file not found:\n%s', baseWaveMcrFile); + end + + data = load(baseWaveMcrFile); + + % Get Hm0 + if isfield(data, 'Hm0') + Hm0 = data.Hm0(:); + elseif isfield(data, 'H') + Hm0 = data.H(:); + elseif isfield(data, 'mcr') && isfield(data.mcr, 'Hm0') + Hm0 = [data.mcr.Hm0].'; + elseif isfield(data, 'mcr') && isfield(data.mcr, 'H') + Hm0 = [data.mcr.H].'; + else + error('Could not find Hm0 or H in base wave-condition file.'); + end + + % Get Tp + if isfield(data, 'Tp') + Tp = data.Tp(:); + elseif isfield(data, 'T') + Tp = data.T(:); + elseif isfield(data, 'mcr') && isfield(data.mcr, 'Tp') + Tp = [data.mcr.Tp].'; + elseif isfield(data, 'mcr') && isfield(data.mcr, 'T') + Tp = [data.mcr.T].'; + else + error('Could not find Tp or T in base wave-condition file.'); + end + + nWave = numel(Hm0); + + % Get Te + if isfield(data, 'Te') + Te = data.Te(:); + elseif isfield(data, 'mcr') && isfield(data.mcr, 'Te') + Te = [data.mcr.Te].'; + else + Te = NaN(nWave, 1); + end + + % Get probability + if isfield(data, 'probability') + probability = data.probability(:); + elseif isfield(data, 'weights') + probability = data.weights(:); + elseif isfield(data, 'mcr') && isfield(data.mcr, 'probability') + probability = [data.mcr.probability].'; + elseif isfield(data, 'mcr') && isfield(data.mcr, 'weights') + probability = [data.mcr.weights].'; + else + probability = ones(nWave, 1) / nWave; + end + + probability = probability / sum(probability); + + % Get condition_id + if isfield(data, 'condition_id') + condition_id = data.condition_id(:); + elseif isfield(data, 'mcr') && isfield(data.mcr, 'condition_id') + condition_id = [data.mcr.condition_id].'; + else + condition_id = (1:nWave).'; + end + + nDamping = numel(dampingValues); + nCases = nWave * nDamping; + + cases = zeros(nCases, 9); + + caseCounter = 0; + + for j = 1:nDamping + for i = 1:nWave + + caseCounter = caseCounter + 1; + + cases(caseCounter, :) = [ ... + caseCounter, ... % 1 caseNum + condition_id(i), ... % 2 condition_id + i, ... % 3 wave_id + j, ... % 4 damping_id + Hm0(i), ... % 5 Hm0 + Tp(i), ... % 6 Tp + Te(i), ... % 7 Te + probability(i), ... % 8 probability + dampingValues(j) ... % 9 damping + ]; + + end + end + + % WEC-Sim-compatible MCR structure + mcr = struct(); + mcr.header = { ... + 'caseNum', ... + 'condition_id', ... + 'wave_id', ... + 'damping_id', ... + 'Hm0', ... + 'Tp', ... + 'Te', ... + 'probability', ... + 'damping' ... + }; + + mcr.cases = cases; + + [runtimeFolder, ~, ~] = fileparts(runtimeMcrFile); + + if ~isempty(runtimeFolder) && ~exist(runtimeFolder, 'dir') + mkdir(runtimeFolder); + end + + save(runtimeMcrFile, 'mcr', 'dampingValues', 'baseWaveMcrFile'); + + fprintf('\nCreated WEC-Sim MCR damping-grid file:\n%s\n', runtimeMcrFile); + fprintf('Wave conditions: %d\n', nWave); + fprintf('Damping values: %d\n', nDamping); + fprintf('Total MCR cases: %d\n\n', nCases); + +end \ No newline at end of file diff --git a/OSWEC_Optimization_Damping/spectrumData.mat b/OSWEC_Optimization_Damping/spectrumData.mat new file mode 100644 index 00000000..4af1086c Binary files /dev/null and b/OSWEC_Optimization_Damping/spectrumData.mat differ diff --git a/OSWEC_Optimization_Damping/userDefinedFunctions.m b/OSWEC_Optimization_Damping/userDefinedFunctions.m new file mode 100644 index 00000000..32ae3d7b --- /dev/null +++ b/OSWEC_Optimization_Damping/userDefinedFunctions.m @@ -0,0 +1,584 @@ +% userDefinedFunctions.m + +%% Detect MCR mode + +isMcrRun = exist('mcr', 'var') && exist('imcr', 'var') && isfield(mcr, 'cases'); + +%% Normal non-MCR WEC-Sim behavior + +if ~isMcrRun + + % Original OSWEC example plots. + % These are intentionally disabled during MCR runs. + + waves.plotElevation(simu.rampTime); + + try + waves.plotSpectrum(); + catch + end + + output.plotForces(1,5) + output.plotResponse(1,5); + output.plotForces(2,1) + + +end + +%% MCR post-processing + +if isMcrRun + + project_folder = fileparts(which('wecSimInputFile.m')); + + resultsFolder = fullfile(project_folder, 'mcr_results'); + + if ~exist(resultsFolder, 'dir') + mkdir(resultsFolder); + end + + caseResultsFile = fullfile(resultsFolder, 'oswec_mcr_case_results.csv'); + caseResultsMat = fullfile(resultsFolder, 'oswec_mcr_case_results.mat'); + summaryFile = fullfile(resultsFolder, 'oswec_damping_optimization_summary.csv'); + summaryMatFile = fullfile(resultsFolder, 'oswec_damping_optimization_summary.mat'); + + % Start fresh on the first MCR case. + if imcr == 1 + + if exist(caseResultsFile, 'file') + delete(caseResultsFile); + end + + if exist(caseResultsMat, 'file') + delete(caseResultsMat); + end + + if exist(summaryFile, 'file') + delete(summaryFile); + end + + if exist(summaryMatFile, 'file') + delete(summaryMatFile); + end + + end + + %% Current MCR case metadata + + caseRow = mcr.cases(imcr, :); + + caseNum = caseRow(1); + condition_id = caseRow(2); + wave_id = caseRow(3); + damping_id = caseRow(4); + Hm0 = caseRow(5); + Tp = caseRow(6); + Te = caseRow(7); + probability = caseRow(8); + damping = caseRow(9); + + %% Extract power and response metrics + + success = true; + error_message = ""; + + try + + [Pmean_W, Pmean_kW, max_pto_power_kW, max_pto_torque_Nm, max_flap_pitch_deg] = ... + oswecExtractMetricsFromBodyPitch(output, damping, simu.rampTime); + + catch ME + + success = false; + error_message = string(ME.message); + + Pmean_W = NaN; + Pmean_kW = NaN; + max_pto_power_kW = NaN; + max_pto_torque_Nm = NaN; + max_flap_pitch_deg = NaN; + + warning('Could not extract metrics for MCR case %d: %s', imcr, ME.message); + + end + + weighted_Pmean_W = probability * Pmean_W; + weighted_Pmean_kW = probability * Pmean_kW; + + %% Store current case result + + thisResult = table( ... + caseNum, ... + condition_id, ... + wave_id, ... + damping_id, ... + Hm0, ... + Te, ... + Tp, ... + probability, ... + damping, ... + success, ... + Pmean_W, ... + Pmean_kW, ... + weighted_Pmean_W, ... + weighted_Pmean_kW, ... + max_pto_power_kW, ... + max_pto_torque_Nm, ... + max_flap_pitch_deg, ... + error_message, ... + 'VariableNames', { ... + 'caseNum', ... + 'condition', ... + 'wave_id', ... + 'damping_id', ... + 'Hm0', ... + 'Te', ... + 'Tp', ... + 'weight', ... + 'damping', ... + 'success', ... + 'Pmean_W', ... + 'Pmean_kW', ... + 'weighted_Pmean_W', ... + 'weighted_Pmean_kW', ... + 'max_pto_power_kW', ... + 'max_pto_torque_Nm', ... + 'max_flap_pitch_deg', ... + 'error_message' ... + } ... + ); + + if exist(caseResultsMat, 'file') + load(caseResultsMat, 'allResults'); + allResults = [allResults; thisResult]; + else + allResults = thisResult; + end + + save(caseResultsMat, 'allResults'); + writetable(allResults, caseResultsFile); + + fprintf('Saved MCR result %d of %d\n', imcr, size(mcr.cases, 1)); + fprintf('Accumulated result rows: %d\n', height(allResults)); + + %% Final MCR case: summarize damping optimization + + if imcr == size(mcr.cases, 1) + + fprintf('\n============================================\n'); + fprintf('All MCR cases complete. Summarizing damping optimization.\n'); + fprintf('============================================\n'); + + optimization_summary = oswecSummarizeDampingResults(allResults); + + writetable(optimization_summary, summaryFile); + + save(summaryMatFile, ... + 'optimization_summary', ... + 'allResults'); + + disp(optimization_summary); + + validRows = optimization_summary.valid_damping == true; + + if any(validRows) + + validSummary = optimization_summary(validRows, :); + + [bestWeightedPowerW, idxBestWeighted] = max(validSummary.weighted_mean_power_W); + bestWeightedDamping = validSummary.damping_Nm_s_per_rad(idxBestWeighted); + + [bestUnweightedPowerW, idxBestUnweighted] = max(validSummary.unweighted_mean_power_W); + bestUnweightedDamping = validSummary.damping_Nm_s_per_rad(idxBestUnweighted); + + bestWeightedPowerkW = bestWeightedPowerW / 1000; + bestUnweightedPowerkW = bestUnweightedPowerW / 1000; + + bestWeightedAEP = bestWeightedPowerkW * 8760; + bestUnweightedAEP = bestUnweightedPowerkW * 8760; + + fprintf('\n--------------------------------------------\n'); + fprintf('Damping optimization complete\n'); + fprintf('--------------------------------------------\n'); + fprintf('Best weighted PTO damping: %.4e N m s/rad\n', bestWeightedDamping); + fprintf('Best weighted mean power: %.3f W\n', bestWeightedPowerW); + fprintf('Best weighted mean power: %.3f kW\n', bestWeightedPowerkW); + fprintf('Best weighted AEP: %.3f kWh/year\n', bestWeightedAEP); + fprintf('--------------------------------------------\n'); + fprintf('Best unweighted PTO damping: %.4e N m s/rad\n', bestUnweightedDamping); + fprintf('Best unweighted mean power: %.3f W\n', bestUnweightedPowerW); + fprintf('Best unweighted mean power: %.3f kW\n', bestUnweightedPowerkW); + fprintf('Best unweighted AEP: %.3f kWh/year\n', bestUnweightedAEP); + fprintf('--------------------------------------------\n'); + + oswecPlotDampingSummary(optimization_summary, resultsFolder); + + else + + warning('No damping values had all successful cases.'); + + end + + fprintf('\nSaved MCR case results:\n%s\n', caseResultsFile); + fprintf('Saved damping optimization summary:\n%s\n', summaryFile); + fprintf('Saved damping optimization MAT:\n%s\n', summaryMatFile); + fprintf('\nAnalysis complete.\n'); + + end + +end + +%% Local helper functions + +function [Pmean_W, Pmean_kW, max_pto_power_kW, max_pto_torque_Nm, max_flap_pitch_deg] = ... + oswecExtractMetricsFromBodyPitch(output, damping, rampTime) +% Extract OSWEC damping-power metrics. +% + + Pmean_W = NaN; + Pmean_kW = NaN; + max_pto_power_kW = NaN; + max_pto_torque_Nm = NaN; + max_flap_pitch_deg = NaN; + + %% Get body 1 output + + if isstruct(output) && isfield(output, 'bodies') + bodyOut = output.bodies(1); + elseif isobject(output) && isprop(output, 'bodies') + bodyOut = output.bodies(1); + else + error('Could not find output.bodies.'); + end + + %% Extract body 1 pitch velocity + + if ~oswecHasMember(bodyOut, 'velocity') + error('Could not find output.bodies(1).velocity.'); + end + + velocityRaw = oswecGetMember(bodyOut, 'velocity'); + + [bodyVelocity, velocityTime] = oswecSignalToArrayAndTime(velocityRaw); + + flapPitchVelocity = oswecExtractDof(bodyVelocity, 5); + + %% Extract body 1 pitch position if available + + flapPitch = []; + + if oswecHasMember(bodyOut, 'position') + + positionRaw = oswecGetMember(bodyOut, 'position'); + + [bodyPosition, positionTime] = oswecSignalToArrayAndTime(positionRaw); + + flapPitch = oswecExtractDof(bodyPosition, 5); + + else + + positionTime = []; + + end + + %% Remove wave ramp period or initial transient + + nSamples = numel(flapPitchVelocity); + + if ~isempty(velocityTime) && numel(velocityTime) == nSamples + + valid = velocityTime >= rampTime; + + else + + % Fallback if no time vector exists: + % ignore first 25 percent of samples. + firstIndex = max(1, floor(0.25 * nSamples)); + valid = false(nSamples, 1); + valid(firstIndex:end) = true; + + end + + flapPitchVelocity = flapPitchVelocity(valid); + + if ~isempty(flapPitch) && numel(flapPitch) == nSamples + flapPitchForMax = flapPitch(valid); + else + flapPitchForMax = flapPitch; + end + + %% Compute PTO power and torque from damping + + ptoPower_W = damping .* flapPitchVelocity.^2; + ptoTorque_Nm = damping .* flapPitchVelocity; + + Pmean_W = mean(ptoPower_W, 'omitnan'); + Pmean_kW = Pmean_W / 1000; + + max_pto_power_kW = max(abs(ptoPower_W), [], 'omitnan') / 1000; + max_pto_torque_Nm = max(abs(ptoTorque_Nm), [], 'omitnan'); + + if ~isempty(flapPitchForMax) + max_flap_pitch_deg = max(abs(flapPitchForMax), [], 'omitnan') * 180 / pi; + end + + %% Diagnostics + + maxVelocity = max(abs(flapPitchVelocity), [], 'omitnan'); + + if isempty(maxVelocity) || isnan(maxVelocity) + error('Flap pitch velocity is empty or NaN after transient removal.'); + end + + if maxVelocity == 0 + warning(['Flap pitch velocity is exactly zero. ', ... + 'Computed absorbed power will be zero. ', ... + 'Check body 1 DOF 5 response and WEC-Sim output fields.']); + end + +end + + +function optimization_summary = oswecSummarizeDampingResults(T) +% Compute optimization metrics by damping. + + damping_values = unique(T.damping); + + weighted_mean_power_W = NaN(size(damping_values)); + unweighted_mean_power_W = NaN(size(damping_values)); + + weighted_mean_power_kW = NaN(size(damping_values)); + unweighted_mean_power_kW = NaN(size(damping_values)); + + weighted_AEP_kWh_per_year = NaN(size(damping_values)); + unweighted_AEP_kWh_per_year = NaN(size(damping_values)); + + valid_damping = false(size(damping_values)); + n_failed = zeros(size(damping_values)); + n_success = zeros(size(damping_values)); + n_total = zeros(size(damping_values)); + + max_flap_pitch_deg_by_damping = NaN(size(damping_values)); + max_pto_torque_Nm_by_damping = NaN(size(damping_values)); + max_pto_power_kW_by_damping = NaN(size(damping_values)); + + for i = 1:length(damping_values) + + damping = damping_values(i); + + rows = T.damping == damping; + T_damping = T(rows, :); + + n_total(i) = height(T_damping); + n_failed(i) = sum(T_damping.success == false); + n_success(i) = sum(T_damping.success == true); + + if n_failed(i) == 0 && n_success(i) > 0 + + weighted_mean_power_W(i) = sum(T_damping.weighted_Pmean_W, 'omitnan'); + unweighted_mean_power_W(i) = mean(T_damping.Pmean_W, 'omitnan'); + + weighted_mean_power_kW(i) = weighted_mean_power_W(i) / 1000; + unweighted_mean_power_kW(i) = unweighted_mean_power_W(i) / 1000; + + weighted_AEP_kWh_per_year(i) = weighted_mean_power_kW(i) * 8760; + unweighted_AEP_kWh_per_year(i) = unweighted_mean_power_kW(i) * 8760; + + valid_damping(i) = true; + + if ismember('max_flap_pitch_deg', T_damping.Properties.VariableNames) + max_flap_pitch_deg_by_damping(i) = max(T_damping.max_flap_pitch_deg, [], 'omitnan'); + end + + if ismember('max_pto_torque_Nm', T_damping.Properties.VariableNames) + max_pto_torque_Nm_by_damping(i) = max(T_damping.max_pto_torque_Nm, [], 'omitnan'); + end + + if ismember('max_pto_power_kW', T_damping.Properties.VariableNames) + max_pto_power_kW_by_damping(i) = max(T_damping.max_pto_power_kW, [], 'omitnan'); + end + + end + + end + + optimization_summary = table( ... + damping_values(:), ... + weighted_mean_power_W(:), ... + weighted_mean_power_kW(:), ... + unweighted_mean_power_W(:), ... + unweighted_mean_power_kW(:), ... + weighted_AEP_kWh_per_year(:), ... + unweighted_AEP_kWh_per_year(:), ... + valid_damping(:), ... + n_success(:), ... + n_failed(:), ... + n_total(:), ... + max_flap_pitch_deg_by_damping(:), ... + max_pto_torque_Nm_by_damping(:), ... + max_pto_power_kW_by_damping(:), ... + 'VariableNames', { ... + 'damping_Nm_s_per_rad', ... + 'weighted_mean_power_W', ... + 'weighted_mean_power_kW', ... + 'unweighted_mean_power_W', ... + 'unweighted_mean_power_kW', ... + 'weighted_AEP_kWh_per_year', ... + 'unweighted_AEP_kWh_per_year', ... + 'valid_damping', ... + 'n_success', ... + 'n_failed', ... + 'n_total', ... + 'max_flap_pitch_deg', ... + 'max_pto_torque_Nm', ... + 'max_pto_power_kW' ... + } ... + ); + + optimization_summary = sortrows(optimization_summary, 'weighted_mean_power_W', 'descend'); + +end + + +function oswecPlotDampingSummary(optimization_summary, resultsFolder) +% Plot final damping optimization summary once. + + valid = optimization_summary.valid_damping == true; + + if ~any(valid) + return + end + + fig = figure; + + semilogx( ... + optimization_summary.damping_Nm_s_per_rad(valid), ... + optimization_summary.weighted_mean_power_kW(valid), ... + 'o', ... + 'LineWidth', 2, ... + 'DisplayName', 'Weighted mean absorbed power'); + + hold on + + semilogx( ... + optimization_summary.damping_Nm_s_per_rad(valid), ... + optimization_summary.unweighted_mean_power_kW(valid), ... + 's', ... + 'LineWidth', 2, ... + 'DisplayName', 'Unweighted mean absorbed power'); + + grid on + xlabel('PTO damping [N m s/rad]') + ylabel('Mean absorbed power [kW]') + title('OSWEC PTO Damping Optimization') + legend('Location', 'best') + + saveas(fig, fullfile(resultsFolder, 'oswec_damping_optimization.png')); + saveas(fig, fullfile(resultsFolder, 'oswec_damping_optimization.fig')); + +end + + +function tf = oswecHasMember(obj, name) +%True if struct field or object property exists. + + if isstruct(obj) + tf = isfield(obj, name); + elseif isobject(obj) + tf = isprop(obj, name); + else + tf = false; + end + +end + + +function value = oswecGetMember(obj, name) +% Get struct field or object property. + + value = obj.(name); + +end + + +function [x, t] = oswecSignalToArrayAndTime(signal) +% Convert common signal formats to numeric array and time. + + t = []; + + if isa(signal, 'timeseries') + t = signal.Time; + x = signal.Data; + x = squeeze(x); + return + end + + if isstruct(signal) + + if isfield(signal, 'time') + t = signal.time; + elseif isfield(signal, 'Time') + t = signal.Time; + end + + if isfield(signal, 'signals') && isfield(signal.signals, 'values') + x = signal.signals.values; + elseif isfield(signal, 'Data') + x = signal.Data; + elseif isfield(signal, 'data') + x = signal.data; + else + x = signal; + end + + else + + x = signal; + + end + + if istable(x) + x = table2array(x); + end + + x = squeeze(x); + +end + + +function dofSignal = oswecExtractDof(x, dof) +%Extract DOF column from WEC-Sim response matrix. +% Handles common layouts: +% N x 6 columns are DOFs +% N x 7 first column is time, columns 2:7 are DOFs +% 6 x N rows are DOFs +% 7 x N first row is time, rows 2:7 are DOFs + + if isvector(x) + + dofSignal = x(:); + return + + end + + [nRows, nCols] = size(x); + + if nCols == 6 + dofSignal = x(:, dof); + + elseif nCols >= 7 + dofSignal = x(:, dof + 1); + + elseif nRows == 6 + dofSignal = x(dof, :).'; + + elseif nRows >= 7 + dofSignal = x(dof + 1, :).'; + + else + error('Could not extract DOF %d from signal with size %d x %d.', dof, nRows, nCols); + end + + dofSignal = dofSignal(:); + +end \ No newline at end of file diff --git a/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_16_clusters.mat b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_16_clusters.mat new file mode 100644 index 00000000..cd920b55 Binary files /dev/null and b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_16_clusters.mat differ diff --git a/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_32_clusters.mat b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_32_clusters.mat new file mode 100644 index 00000000..52a66bd4 Binary files /dev/null and b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_32_clusters.mat differ diff --git a/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_4_clusters.mat b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_4_clusters.mat new file mode 100644 index 00000000..875df9d6 Binary files /dev/null and b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_4_clusters.mat differ diff --git a/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_64_clusters.mat b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_64_clusters.mat new file mode 100644 index 00000000..d2391e0e Binary files /dev/null and b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_64_clusters.mat differ diff --git a/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_8_clusters.mat b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_8_clusters.mat new file mode 100644 index 00000000..696000ed Binary files /dev/null and b/OSWEC_Optimization_Damping/wave_conditions/wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_8_clusters.mat differ diff --git a/OSWEC_Optimization_Damping/wecSimInputFile.m b/OSWEC_Optimization_Damping/wecSimInputFile.m new file mode 100644 index 00000000..69555df4 --- /dev/null +++ b/OSWEC_Optimization_Damping/wecSimInputFile.m @@ -0,0 +1,128 @@ +% wecSimInputFile.m +% OSWEC input file for WEC-Sim MCR damping sweep. +% 1. Set dampingValues below. +% 2. Run: wecSimMCR + +%% MCR Damping-Sweep Setup + +project_folder = fileparts(mfilename('fullpath')); + +% Base wave-condition file generated by Python/MHKiT. +% This can be either your array-style .mat or your wave-condition mcr .mat. +baseWaveMcrFile = fullfile(project_folder, 'wave_conditions', ... + 'wave_conditions_buoy_46050_2020_2021_2022_2023_2024_2025_16_clusters.mat'); + +% Runtime WEC-Sim-compatible MCR file. +runtimeMcrFile = fullfile(project_folder, 'wave_conditions', ... + 'oswec_runtime_damping_grid_mcr.mat'); + +% PTO damping values for the grid search. +% Change this line to change the damping sweep. +dampingValues = [1e4 1e5 1e6 1e7 1e8].'; + +% Build the runtime MCR file on the initial setup evaluation. +% During the actual MCR loop, variable mcr should already exist. +if ~exist('mcr', 'var') || ~isfield(mcr, 'cases') + + runtimeMcrFile = oswecBuildDampingMCR( ... + baseWaveMcrFile, ... + dampingValues, ... + runtimeMcrFile); + + data = load(runtimeMcrFile, 'mcr'); + mcr = data.mcr; + +end + +% Determine current MCR case index. +% In WEC-Sim MCR, the loop variable is commonly imcr. +if exist('imcr', 'var') + caseIndex = imcr; +elseif isfield(mcr, 'caseNum') + caseIndex = mcr.caseNum; +else + caseIndex = 1; +end + +% Pull current case row from mcr.cases. +caseRow = mcr.cases(caseIndex, :); + +% Column mapping: +caseNum = caseRow(1); +conditionID = caseRow(2); +waveID = caseRow(3); +dampingID = caseRow(4); +caseHm0 = caseRow(5); +caseTp = caseRow(6); +caseTe = caseRow(7); +caseProb = caseRow(8); +caseDamping = caseRow(9); + +%% Simulation Data + +simu = simulationClass(); + +% Use state-space radiation formulation +simu.stateSpace = 0; + +simu.simMechanicsFile = 'OSWEC.slx'; +simu.mode = 'normal'; +simu.explorer = 'off'; + +% WEC-Sim-compatible MCR file +simu.mcrMatFile = runtimeMcrFile; + +simu.startTime = 0; +simu.rampTime = 50; +simu.endTime = 400; + +simu.solver = 'ode45'; +simu.dt = 0.01; + +simu.cicEndTime = 30; + +%% Wave Information + +waves = waveClass('irregular'); + +waves.height = caseHm0; % Significant wave height Hm0 [m] +waves.period = caseTp; % Peak period Tp [s] +waves.spectrumType = 'JS'; % JONSWAP spectrum + +%% Body Data + +body(1) = bodyClass('../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); +body(1).geometryFile = '../_Common_Input_Files/OSWEC/geometry/flap.stl'; +body(1).mass = 127000; +body(1).inertia = [1.85e6 1.85e6 1.85e6]; + +body(2) = bodyClass('../_Common_Input_Files/OSWEC/hydroData/oswec.h5'); +body(2).geometryFile = '../_Common_Input_Files/OSWEC/geometry/base.stl'; +body(2).mass = 999; +body(2).inertia = [999 999 999]; + +%% PTO and Constraint Parameters + +constraint(1) = constraintClass('Constraint1'); +constraint(1).location = [0 0 -10]; + +pto(1) = ptoClass('PTO1'); +pto(1).stiffness = 0; +pto(1).damping = caseDamping; +pto(1).location = [0 0 -8.9]; + +%% Console Output + +fprintf('\n--------------------------------------------\n'); +fprintf('OSWEC MCR damping-sweep case\n'); +fprintf('MCR case number: %d\n', caseNum); +fprintf('Wave condition ID: %d\n', conditionID); +fprintf('Wave ID: %d\n', waveID); +fprintf('Damping ID: %d\n', dampingID); +fprintf('Hm0 = %.3f m\n', caseHm0); +fprintf('Tp = %.3f s\n', caseTp); +fprintf('Te = %.3f s\n', caseTe); +fprintf('Probability = %.5f\n', caseProb); +fprintf('PTO damping = %.3e\n', caseDamping); +fprintf('MCR file: %s\n', runtimeMcrFile); +fprintf('--------------------------------------------\n\n'); \ No newline at end of file