MHKiT Loads Module

The following example will help familiarize you with some of the functions in the MHKiT loads module that you can use to assist you in your loads analysis.

Start by importing the necessary python packages and MHKiT module.

[1]:
import pandas as pd
import numpy as np
from mhkit import utils
from mhkit import loads
import matplotlib.pyplot as plt

Import Loads Data

This example data comes from a land based wind turbine, since there is limited availablality of loads data for MHK devices. The example uses a subset of data from a database of 331 files, each containing 10 minutes of data sampled at 50Hz.

As a start, lets look at the data for one of these files to figure out what formatting we need to apply. We utilize pandas to read in the csv and peek into what the data looks like.

[2]:
loads_data_file = "./data/loads/data_loads_example.csv"

# Import csv data file
raw_loads_data = pd.read_csv(loads_data_file)
raw_loads_data.head()
[2]:
Timestamp Time uWind_80m WD_ModActive WD_Nacelle WD_NacelleMod LSSDW_Tq LSSDW_My LSSDW_Mz TTTq TT_ForeAft TT_SideSide TB_ForeAft TB_SideSide BL3_FlapMom BL3_EdgeMom BL1_FlapMom BL1_EdgeMom ActivePower yawoffset
0 42795.061586 0.00 3.226754 1.0 157.302829 157.279582 -41.380694 -234.487436 -6.207381 -70.130726 -936.247028 -12.605151 -330.410779 1024.816867 470.774738 -165.541786 33.427748 -59.452360 -5.279680 0.023247
1 42795.061586 0.02 3.221099 1.0 157.302829 157.279582 -38.614459 -233.715870 -8.886200 -66.916338 -942.675906 -24.350452 -315.445562 873.214212 469.244736 -163.588005 32.697822 -62.300637 -5.671178 0.023247
2 42795.061586 0.04 3.223492 1.0 157.302829 157.279582 -39.717967 -234.341966 -7.970862 -67.860011 -922.971018 -22.485796 -292.252115 876.299461 468.736474 -166.018111 35.495810 -61.733604 -5.551847 0.023247
3 42795.061586 0.06 3.223274 1.0 157.302829 157.279582 -41.415143 -235.645598 -10.451794 -72.371951 -939.515265 -33.030892 -274.812769 763.833828 467.373452 -164.645639 37.952455 -64.390050 -4.626557 0.023247
4 42795.061587 0.08 3.223927 1.0 157.302829 157.279582 -38.614459 -234.755991 -8.648988 -76.530014 -924.771486 -29.228398 -310.213400 704.537757 466.318754 -161.233863 37.430668 -65.974766 -4.708621 0.023247

Format Loads Data with datetime index

To use MHKiT it is important to index your DataFrame by datetime. The example loads data has two references to time, but neither are in the right format. The 'Timestamp' column is what will give us the datetime index that we need, but we first need to convert it from microsoft excel format to pd.Datetime.

[3]:
# Use the datetime conversion from the utils module
datetime = utils.excel_to_datetime(raw_loads_data["Timestamp"])

# Replace the 'Timestamp' column with our newly formatted datetime
raw_loads_data["Timestamp"] = datetime

# Set this as our index for our DataFrame
loads_data = raw_loads_data.set_index("Timestamp")

# Remove the 'time' column since it will not be used
loads_data.drop(columns="Time", inplace=True)
loads_data.head()
[3]:
uWind_80m WD_ModActive WD_Nacelle WD_NacelleMod LSSDW_Tq LSSDW_My LSSDW_Mz TTTq TT_ForeAft TT_SideSide TB_ForeAft TB_SideSide BL3_FlapMom BL3_EdgeMom BL1_FlapMom BL1_EdgeMom ActivePower yawoffset
Timestamp
2017-03-01 01:28:40.999987200 3.226754 1.0 157.302829 157.279582 -41.380694 -234.487436 -6.207381 -70.130726 -936.247028 -12.605151 -330.410779 1024.816867 470.774738 -165.541786 33.427748 -59.452360 -5.279680 0.023247
2017-03-01 01:28:41.020032000 3.221099 1.0 157.302829 157.279582 -38.614459 -233.715870 -8.886200 -66.916338 -942.675906 -24.350452 -315.445562 873.214212 469.244736 -163.588005 32.697822 -62.300637 -5.671178 0.023247
2017-03-01 01:28:41.039990400 3.223492 1.0 157.302829 157.279582 -39.717967 -234.341966 -7.970862 -67.860011 -922.971018 -22.485796 -292.252115 876.299461 468.736474 -166.018111 35.495810 -61.733604 -5.551847 0.023247
2017-03-01 01:28:41.060035200 3.223274 1.0 157.302829 157.279582 -41.415143 -235.645598 -10.451794 -72.371951 -939.515265 -33.030892 -274.812769 763.833828 467.373452 -164.645639 37.952455 -64.390050 -4.626557 0.023247
2017-03-01 01:28:41.079993600 3.223927 1.0 157.302829 157.279582 -38.614459 -234.755991 -8.648988 -76.530014 -924.771486 -29.228398 -310.213400 704.537757 466.318754 -161.233863 37.430668 -65.974766 -4.708621 0.023247

Loads Analysis

Now that we have our loads data in the correct format for MHKiT, we do some analysis.

Damage Equivalent Loads

Lets say that we wanted to investigate fatigue. We can do this by calculating short-term damage equivalent loads (DELs). In this instance, we calculate the DELs on the tower base moment 'TB_ForeAft', and on blade 1 root flap moment 'BL1_FlapMom'. Our tower is steel while our blade is composite so they will have different material slopes. We will run the function damage_equivalent_load on each data signal.

We call our function and apply the default inputs of using at least 100 bins for the load ranges and we let data_length=600 seconds so that we get an equivalent 1Hz DEL for our 10 minute file.

[4]:
# Calculate the damage equivalent load for blade 1 root momement and tower base moment
DEL_tower = loads.general.damage_equivalent_load(
    loads_data["TB_ForeAft"], 4, bin_num=100, data_length=600
)
DEL_blade = loads.general.damage_equivalent_load(
    loads_data["BL1_FlapMom"], 10, bin_num=100, data_length=600
)
print("DEL TB_ForeAft: " + str(DEL_tower))
print("DEL BL1_FlapMom: " + str(DEL_blade))
DEL TB_ForeAft: 3912.6390862773
DEL BL1_FlapMom: 1435.8222478714554

Calculate Statistics

Another important part of loads analysis is looking at statistics. Here, we use another function to help us calculate the mean, max, min, and std for this 10 minute file. Per standards, a valid statistical window has to be consecutive in time with the correct number of datapoints. If this 10 minute file did meet this criteria, then no stats would be generated and a warning message would appear.

NOTE: Sometimes individual files may contain enough data for multiple statistical windows. This function can still handle this scenario as long as the correct inputs are specified.

[5]:
# Calculate the means, maxs, mins, and stdevs for all data signals in the loads data file
means, maxs, mins, stdevs = utils.get_statistics(loads_data, 50, period=600)

# Display the results, indexed by the first timestamp of the corresponding statistical window
means
[5]:
uWind_80m WD_ModActive WD_Nacelle WD_NacelleMod LSSDW_Tq LSSDW_My LSSDW_Mz TTTq TT_ForeAft TT_SideSide TB_ForeAft TB_SideSide BL3_FlapMom BL3_EdgeMom BL1_FlapMom BL1_EdgeMom ActivePower yawoffset
2017-03-01 01:28:41 7.773325 1.0 178.612256 178.602595 127.244191 -252.23813 3.50322 7.032573 -846.663367 271.446574 3785.034515 7.199176 -494.858287 266.790368 -452.652744 21.259999 234.578289 0.009661

At this point, it would be nice to start visualizing some of this data. In order to do this, we need to calculate the stats and DELs for all the files in our database. In this case, it would be done through a loop that imports each file and applies all the functions we just saw. At the end of the loop, we store the result by appending to a list (this is computationally more efficient than concatenating dataframes). Finally, we can convert our lists to dataframes so that its easier to play with the data. To speed things up, this was already done so we just need to import the resulting dataframes.

As a reference, an example code block is shown of how to create a loop that calculates all the means for each file which are then stored into a dataframe.

# pre-allocate lists for storage
means = []
time = []

# start loop
for f in os.listdir(pathOut):
    if f.endswith('.csv'):
        # import csv file
        raw_loads_data = pd.read_csv(pathOut+'/'+f)
        # replace the timestamp column with formatted datetime
        datetime = utils.excel_to_datetime(raw_loads_data['Timestamp'])
        raw_loads_data['Timestamp'] = datetime
        # set this as our index for our dataframe
        loads_data = raw_loads_data.set_index('Timestamp')
        # remove the "time" column as its unnecessary
        loads_data.drop(columns='Time',inplace=True)
        # get stats
        fmean, fmax, fmin, fstd = utils.get_statistics(loads_data,freq=50,period=600)
        means.append(fmean.values.tolist())
        time.append(fmean.index.values)

# convert lists into a dataframe
loads_means = pd.DataFrame(np.squeeze(means),columns=loads_data.columns.values,index=time)
[6]:
# Load DataFrames containing load statistics
means = pd.read_csv("./data/loads/data_loads_means.csv")
maxs = pd.read_csv("./data/loads/data_loads_maxs.csv")
mins = pd.read_csv("./data/loads/data_loads_mins.csv")
std = pd.read_csv("./data/loads/data_loads_std.csv")

means.head()
[6]:
uWind_80m WD_ModActive WD_Nacelle WD_NacelleMod LSSDW_Tq LSSDW_My LSSDW_Mz TTTq TT_ForeAft TT_SideSide TB_ForeAft TB_SideSide BL3_FlapMom BL3_EdgeMom BL1_FlapMom BL1_EdgeMom ActivePower yawoffset
0 7.773325 1.0 178.612256 178.602595 127.244191 -252.238130 3.503220 7.032573 -846.663367 271.446574 3785.034515 7.199176 -494.858288 266.790368 -452.652743 21.259999 234.578289 0.009661
1 4.294855 1.0 171.095503 171.104399 6.705063 -242.954279 8.781801 13.574055 -1005.504041 210.432881 504.631955 -63.459058 -123.664013 179.906995 -95.184346 -44.686204 32.156167 -0.008896
2 5.210606 1.0 168.688106 168.680758 51.782698 -224.455118 15.310279 -53.986054 -926.132830 138.439249 1992.287971 -42.414665 -314.992337 228.895991 -270.323145 -5.000843 87.337237 0.007348
3 14.210652 1.0 182.695493 182.688322 386.484523 -218.151791 13.457946 -39.223265 -568.024026 515.128540 7075.518787 649.320139 -689.885628 301.399047 -718.343755 77.781426 692.061262 0.007171
4 10.558234 1.0 182.443087 182.426475 561.122866 -228.167278 -30.095950 -51.686133 -347.327549 621.974202 10992.154570 598.674976 -1089.599789 374.679600 -1151.084828 138.214755 997.975514 0.016613

Plot Statistics

Now that we have the load statistics, lets display the data as scatter plot. Using the plot_statistics function, we can quickly create a standard scatter plot showing how load variables trend with wind speed. Using this we can quickly identify expected trends and track down outliers.

[7]:
loads.graphics.plot_statistics(
    means["uWind_80m"],
    means["BL1_FlapMom"],
    maxs["BL1_FlapMom"],
    mins["BL1_FlapMom"],
    y_stdev=std["BL1_FlapMom"],
    xlabel="Wind Speed [m/s]",
    ylabel="Blade Flap Moment [kNm]",
    title="Blade Flap Moment Load Statistics",
)

loads.graphics.plot_statistics(
    means["uWind_80m"],
    means["TB_ForeAft"],
    maxs["TB_ForeAft"],
    mins["TB_ForeAft"],
    y_stdev=std["TB_ForeAft"],
    xlabel="Wind Speed [m/s]",
    ylabel="Tower Base Moment [kNm]",
    title="Tower Base Moment Load Statistics",
)
_images/loads_example_13_0.png
_images/loads_example_13_1.png
[7]:
<AxesSubplot:title={'center':'Tower Base Moment Load Statistics'}>

Another common step is to bin the statistical data. This can easily be done with the bin_stats function from the loads module shown below. A warning message will show if there are any bins that were not filled.

[8]:
# Create array containing wind speeds to use as bin edges
bin_edges = np.arange(3, 26, 1)
bin_against = means["uWind_80m"]

# Apply function for means, maxs, and mins
[bin_means, bin_means_std] = loads.general.bin_statistics(means, bin_against, bin_edges)
[bin_maxs, bin_maxs_std] = loads.general.bin_statistics(maxs, bin_against, bin_edges)
[bin_mins, bin_mins_std] = loads.general.bin_statistics(mins, bin_against, bin_edges)

bin_means
Warning: some bins may be empty!
Warning: some bins may be empty!
Warning: some bins may be empty!
[8]:
uWind_80m WD_ModActive WD_Nacelle WD_NacelleMod LSSDW_Tq LSSDW_My LSSDW_Mz TTTq TT_ForeAft TT_SideSide TB_ForeAft TB_SideSide BL3_FlapMom BL3_EdgeMom BL1_FlapMom BL1_EdgeMom ActivePower yawoffset
0 3.582788 1.000000 179.116583 179.215605 29.890756 -252.235192 -2.931479 -26.271654 -932.953738 126.807011 994.638671 -102.050137 -206.591291 129.390923 -238.085263 -14.874877 21.370858 -0.099022
1 4.498995 1.000000 178.968041 178.963826 77.559257 -248.081262 -3.389182 -51.890684 -856.012087 154.470046 2235.672810 55.864708 -327.707803 162.318746 -354.474414 -7.655706 76.033527 0.004216
2 5.525786 0.980392 180.016225 179.996997 171.459516 -230.697510 -2.588070 -50.333653 -761.047719 213.700239 3822.643018 137.424215 -495.457116 191.503291 -599.658658 1436.266991 167.928433 0.019228
3 6.520641 0.970588 180.835469 180.788600 219.138550 -237.979417 -5.688640 -42.827109 -647.310249 302.591017 6122.425105 192.772384 -708.321569 237.072176 -777.545681 599.754365 317.629187 0.046869
4 7.534458 1.000000 181.373937 181.474371 285.072650 -240.838000 -8.602511 -56.222891 -571.093365 376.755167 7551.899048 296.723677 -861.049383 252.077699 -874.336546 88.957621 454.388684 -0.100434
5 8.483451 0.966667 180.006753 180.002341 405.504595 -233.796744 -7.596248 -32.233513 -430.681550 502.070942 9938.523640 335.667789 -1073.556266 313.335072 -1167.966930 1302.506123 679.886822 0.004413
6 9.635533 0.880000 178.894596 178.885289 581.709245 -218.693601 -8.801577 -72.657764 -330.525678 668.243228 11477.966048 588.384954 -1209.842191 342.343342 -1526.574097 4686.436750 942.635185 0.009307
7 10.552790 0.892857 179.473587 179.487138 653.097969 -224.998487 -5.756869 -62.714160 -283.620699 750.186980 12002.102796 740.834502 -1234.723513 354.710593 -1380.119025 2065.459821 1101.930764 -0.013551
8 11.431282 0.814807 179.157452 179.173400 738.418861 -209.685631 -1.320960 -22.419344 -295.116952 874.826602 12237.843754 668.581523 -1248.802900 362.942457 -1291.630102 173.010056 1251.222715 -0.015947
9 12.447532 0.888889 178.668483 178.658239 791.696059 -188.741923 -2.588168 -29.769987 -291.021580 902.150404 11648.909089 743.454123 -1173.606189 328.218911 -1262.103186 738.232546 1316.246001 0.010243
10 13.406650 0.833333 180.438968 180.452124 817.388823 -198.989294 -7.197755 -24.899471 -331.564307 992.711011 10733.428547 791.711119 -1074.877395 300.227189 -1114.404530 121.228715 1384.312354 -0.013155
11 14.348479 0.542676 179.354862 179.346215 875.198293 -59.267478 -1.437998 4.153358 -374.601981 907.552054 9896.827494 677.680004 -971.840066 285.582399 -1036.798595 102.258425 1350.838415 0.008648
12 15.579021 0.833333 178.013520 178.001359 880.603457 -164.536103 -1.851375 -26.848873 -330.960964 1059.024803 9234.857465 926.368719 -899.553762 227.524012 -923.851217 74.708350 1491.638595 0.012160
13 16.400735 1.000000 178.314956 178.311169 866.227371 -185.673240 -12.999453 -78.397677 -287.228505 1071.363839 8983.182063 984.403585 -860.376373 200.755492 -864.637118 53.896186 1499.435045 0.003787
14 17.602406 1.000000 179.171230 179.157532 886.425220 -137.146930 -8.231045 -3.156350 -383.416299 1063.362197 8563.213120 873.449181 -762.282362 225.211156 -810.106716 48.225543 1476.244726 0.013698
15 18.982808 1.000000 179.199942 179.421965 863.882390 -158.281553 -23.757701 -74.334173 -270.882688 1056.399578 8206.587839 1349.955119 -699.125685 205.630765 -703.231687 14.362069 1514.168487 -0.222023
16 19.519921 1.000000 178.361870 178.356581 881.739431 -169.696939 1.018761 -61.347061 -284.548611 1113.928290 7807.792946 981.873207 -695.129807 151.699529 -709.755338 20.152316 1514.549099 0.005289
17 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
18 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
19 22.062534 1.000000 178.524018 178.507149 854.648338 -142.928299 -43.759922 -111.305482 -262.456906 1067.274340 7566.500039 1369.827084 -572.629544 186.356598 -575.251752 4.025634 1500.616164 0.016870
20 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
21 24.759699 1.000000 179.912515 179.904301 334.771556 -206.441986 -11.694631 -32.524912 -679.532690 492.819279 3419.356234 752.087159 -226.330646 137.392161 -187.876792 -33.225278 615.054202 0.008214

Now lets make some more plots with the binned data. Here we use the binned data and corresponding standard deviations as inputs to the plot_bin_statistics function.

[9]:
# Specify center of each wind speed bin, and signal name for analysis
bin_centers = np.arange(3.5, 25.5, step=1)
signal_name = "TB_ForeAft"

# Specify inputs to be used in plotting
bin_mean = bin_means[signal_name]
bin_max = bin_maxs[signal_name]
bin_min = bin_mins[signal_name]
bin_mean_std = bin_means_std[signal_name]
bin_max_std = bin_maxs_std[signal_name]
bin_min_std = bin_mins_std[signal_name]

# Plot binned statistics
loads.graphics.plot_bin_statistics(
    bin_centers,
    bin_mean,
    bin_max,
    bin_min,
    bin_mean_std,
    bin_max_std,
    bin_min_std,
    xlabel="Wind Speed [m/s]",
    ylabel=signal_name,
    title="Binned Statistics",
)
_images/loads_example_17_0.png
[9]:
<AxesSubplot:title={'center':'Binned Statistics'}>