NREL5W_TCN_DDModel_UniWind_rev04.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"collapsed_sections":[],"machine_shape":"hm","mount_file_id":"1bwqApjqzSFmPZjtNiAwj3_i5ns41hVBS","authorship_tag":"ABX9TyPlGwBves2Txcy0zmo6K626"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"},"accelerator":"GPU","gpuClass":"standard"},"cells":[{"cell_type":"markdown","metadata":{"id":"03085429-d9ff-4746-a9a1-5240e9e945d3"},"source":["# NREL 5MW Temporal Conv Net (TCN) Surrogate model\n","\n","This notebook is the initial try to build a TCN surrogate model on NREL 5MW aero-elastic simulations.\n","\n","\n","*   The wind is uniform over the rotor, the input to the model is HH wind time seris (U)\n","*   Normalization using MaxMin scaling\n","\n","[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/googlecolab/colabtools/blob/master/notebooks/colab-github-demo.ipynb)\n","\n"]},{"cell_type":"markdown","source":["# Mounting Google Drive\n","All the simulations are save in CSV files and saved on the google drive."],"metadata":{"id":"WXwFoJQHXu6Q"}},{"cell_type":"markdown","metadata":{"id":"e61e99b6-00f9-4653-8959-be23050c26d2"},"source":["## Installing and Loading packages\n","Some packages are required for this work. Mainly, the TCN package developed by  Philippe Rémy (https://github.com/philipperemy/keras-tcn)."]},{"cell_type":"code","source":["from google.colab import drive\n","import os\n","drive.mount('/content/drive/')\n","\n"],"metadata":{"id":"73TRW8XMYBCg","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1665180703514,"user_tz":420,"elapsed":1035,"user":{"displayName":"Rad Haghi","userId":"05221314968911428197"}},"outputId":"5814cbe4-34b2-4624-95bb-7e3d10f200bb"},"execution_count":5,"outputs":[{"output_type":"stream","name":"stdout","text":["Drive already mounted at /content/drive/; to attempt to forcibly remount, call drive.mount(\"/content/drive/\", force_remount=True).\n"]}]},{"cell_type":"markdown","source":["## Change the current working directory\n","To make life easier, the current directory changed to the folder."],"metadata":{"id":"1F0clXDvanPL"}},{"cell_type":"code","source":["os.chdir('/content/drive/MyDrive/Colab Notebooks/NREL5MW_TCN_DD_UniWind')\n","print(os.getcwd())"],"metadata":{"id":"dbcdkPe3atG-","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1665180706097,"user_tz":420,"elapsed":163,"user":{"displayName":"Rad Haghi","userId":"05221314968911428197"}},"outputId":"6baa1132-0c8f-4dc8-cbf2-0c8c9b10bfba"},"execution_count":6,"outputs":[{"output_type":"stream","name":"stdout","text":["/content/drive/MyDrive/Colab Notebooks/NREL5MW_TCN_DD_UniWind\n"]}]},{"cell_type":"markdown","metadata":{"id":"9d16c926-ef68-4071-b93a-fd5658b6e49a"},"source":["KERAS-TCN installation"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"e8d6c75f-f031-4a55-b88e-836016284b24"},"outputs":[],"source":["!pip3 install keras-tcn"]},{"cell_type":"markdown","source":["I like to have a latex compiler for the figures. Here is the package."],"metadata":{"id":"S-RQeXSt62Il"}},{"cell_type":"code","source":["!sudo apt install cm-super dvipng texlive-latex-extra texlive-latex-recommended\n"],"metadata":{"id":"4AoLFkeY1Zix"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"c08d9a88-9d6a-42e0-aa00-79872d7e769d"},"source":["Other required packges are called in this cell. Also the fonts are set for the figures."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"9bca1312-04a1-4ec7-9ab9-c1e0ee8f1cbe"},"outputs":[],"source":["import tensorflow as tf\n","import os\n","import numpy\n","import time\n","from utilities import select_gpus, plot_results # utilities.py: Contains a few miscellaneous functions \n","#from tcnae import TCNAE # tcnae.py: Specification of the TCN-AE model\n","#import data # data.py: Allows to generate anomalous Mackey-Glass (MG) time series\n","import pandas as pd\n","import numpy as np\n","from sklearn.preprocessing import MinMaxScaler\n","from sklearn.preprocessing import StandardScaler\n","import matplotlib.pyplot as plt\n","import h5py\n","\n","\n","import time as tm\n","import tensorflow\n","from tcn import TCN, tcn_full_summary\n","\n","\n","from tensorflow.keras.layers import Dense, Activation, Dropout\n","from tensorflow.keras.layers import Input\n","from tensorflow.keras.layers import Conv1D\n","from tensorflow.keras.layers import UpSampling1D\n","from tensorflow.keras.layers import AveragePooling1D\n","from tensorflow.keras.layers import MaxPooling1D\n","from tensorflow.keras import optimizers\n","from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint\n","from tensorflow.keras.models import Model\n","from tensorflow.keras import Sequential\n","\n","import random\n","\n","\n","\n","plt.rc('text', usetex=True)\n","plt.rc('font', family='serif')\n","plt.rc('font', size = 30)\n","plt.rc('xtick', labelsize = 25)\n","plt.rc('ytick', labelsize = 25)\n","plt.rc('legend',fontsize=25) # using a size in points\n"]},{"cell_type":"markdown","source":["## Test GPU existance and employe one"],"metadata":{"id":"rVXHycNxcrqq"}},{"cell_type":"code","source":["device_name = tf.test.gpu_device_name()\n","if device_name != '/device:GPU:0':\n","  raise SystemError('GPU device not found')\n","print('Found GPU at: {}'.format(device_name))"],"metadata":{"id":"DOqJZ4fMcxVO"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"5e5c7ca9-053c-426c-9af4-c5edbc7c5160"},"source":["\n","\n","# Data import from CSV saved before\n","\n","Read the data channels and save them in a `dictionary` format."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"e4b24422-36c1-4b47-a966-7478aa3280f2"},"outputs":[],"source":["# Reading the wind speeds\n","raw_data = {}\n","raw_data['wspdX'] = pd.read_csv('data_fullDLC12/WindSpeedX.csv').transpose()\n","raw_data['MyB'] = pd.read_csv('data_fullDLC12/RootMyb.csv').transpose()\n","raw_data['MxB'] = pd.read_csv('data_fullDLC12/RootMxb.csv').transpose()\n","raw_data['MxTB'] = pd.read_csv('data_fullDLC12/TwrBsMy.csv').transpose()\n","raw_data['MyTB'] = pd.read_csv('data_fullDLC12/TwrBsMy.csv').transpose()\n","raw_data['FxTT'] = pd.read_csv('data_fullDLC12/YawBrFxp.csv').transpose()\n","raw_data['FyTT'] = pd.read_csv('data_fullDLC12/YawBrFyp.csv').transpose()\n","raw_data['MxTT'] = pd.read_csv('data_fullDLC12/YawBrMxp.csv').transpose()\n","raw_data['MyTT'] = pd.read_csv('data_fullDLC12/YawBrMyp.csv').transpose()\n","\n","\n","\n","\n"]},{"cell_type":"markdown","source":["# Nested dictionaries\n","Building the nested dictionries from the importated data, and normalizing the input-output data based on the MaxMin scaling"],"metadata":{"id":"5GRc78tWl9Is"}},{"cell_type":"code","source":["chls = ['wspdX','MxB','MyB','FxTT','FyTT','MxTT','MyTT','MxTB','MyTB']\n","mps = ['03mps','04mps','05mps','06mps','07mps','08mps','09mps','10mps','11mps','12mps','13mps','14mps','15mps','16mps','17mps','18mps','19mps','20mps','21mps','22mps','23mps','24mps','25mps']\n","\n","data_sets={}\n","data_sets_norm={}\n","\n","for ch in chls:\n","  inx_lst = raw_data[ch].index.to_list()\n","  data_sets[ch] = {}\n","  data_sets_norm[ch] ={}\n","  for u in mps:\n","    indices = [i for i, s in enumerate(inx_lst) if u in s]\n","    data_sets[ch][u] = raw_data[ch].iloc[indices]\n","    data_sets_norm[ch][u] =data_sets[ch][u].subtract(data_sets[ch][u].min().min(),axis=0).div((data_sets[ch][u].max().max()-data_sets[ch][u].min().min()),axis=0)\n","    del indices\n","\n","\n","\n"],"metadata":{"id":"QmeFzf3Dx7s8"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"93e47f6d-0a3d-4063-be7b-791c9c72af7e"},"source":["#### Output normalization visual check"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"ea7bace8-20bb-4496-96f2-66e57d3c9b29"},"outputs":[],"source":["u = '22mps'\n","idx = 12 #This is an arbitary number to present one of the seeds as an example. It can be anything between 0-90.\n","\n","fig,ax = plt.subplots(4,1,sharex=True)\n","ax[0].plot(data_sets_norm['wspdX'][u].iloc[idx])\n","ax[1].plot(data_sets_norm['MxB'][u].iloc[idx])\n","ax[2].plot(data_sets_norm['MyB'][u].iloc[idx])\n","ax[3].plot(data_sets_norm['MxTB'][u].iloc[idx])\n"]},{"cell_type":"markdown","metadata":{"id":"4657d7ea-857a-4e52-9d31-4770b1177f47"},"source":["# TCN-AE trial for the wind input and loads output\n","\n","What I am going to try here is training the TCN-AE network on the different input output instead of the same input-output. I want to see if things are going to work. "]},{"cell_type":"code","execution_count":null,"metadata":{"id":"9bc53815-f347-4095-a16d-cfd6067f7222"},"outputs":[],"source":["import tcnae\n","import imp\n","imp.reload(tcnae)\n","\n","# This is for the naming convention of the saved figures\n","timestr = tm.strftime('%Y%m%d')\n","\n","# Time for x axis\n","time = np.linspace(0,600,12000)\n","\n","# For every wind speed we have 90 seeds\n","seeds = np.arange(0,90)\n","\n","# The x-axis for the loss plot\n","epoch_xaxis=np.arange(1,301)\n","\n","# Num of seeds used for the training\n","TrnSzSd = 40\n","\n","# The wind speeds the training is going to be done\n","U=['24mps']\n","\n","# The channels/sensors the training is going to be done\n","CH=['MyTB']\n","\n","# Dict to save the history of training\n","hist_data_set = {}\n","\n","#Dict to save the the predictions values\n","ypred={}\n","\n","for u in U:\n","  hist_data_set[u] = {}\n","  ypred[u]={}\n","  for ch in CH:\n","    print(u,ch)\n","    # Choose seeds out of 90 seeds for training\n","    rnd_seed = np.random.random_integers(1, high=90, size=TrnSzSd)\n","    # Find the remaninag seeds and use them later for prediction and choose \n","    # one random\n","    rem_seeds = seeds[~np.in1d(np.arange(seeds.size),rnd_seed)]\n","    SeedNo=random.choice(rem_seeds)\n","    # Dedicate input and output to x and y for training\n","    x = data_sets_norm['wspdX'][u].iloc[rnd_seed,1:]\n","    y = data_sets_norm[ch][u].iloc[rnd_seed,1:]\n","    # Build the model based on the setting in the 'tcnae.py' file\n","    tcn_ae = tcnae.TCNAE()\n","    tcn_ae.build_model()\n","    # Fit the model to the traininag data\n","    history = tcn_ae.model.fit(x, y, \n","                            batch_size=12, \n","                            validation_split=0.1, \n","                            shuffle=True,\n","                            verbose=1,\n","                            epochs = 300)\n","    # Keep thr history\n","    hist_data_set[u][ch] = history\n","    # Prediction with the random seed\n","    y_pred = tcn_ae.model.predict(data_sets_norm['wspdX'][u].iloc[SeedNo:SeedNo+1,1:])\n","    # Saving the prediction\n","    ypred[u][ch] = y_pred\n","    # Plotting\n","    fig = plt.figure(figsize=(20,14))\n","    gs = fig.add_gridspec(2,2)\n","    ax1 = fig.add_subplot(gs[0, :])\n","    ax2 = fig.add_subplot(gs[1, 0])\n","    ax3 = fig.add_subplot(gs[1, 1])\n","    ax1.semilogy(epoch_xaxis,history.history['loss'],linewidth=2,label='Training Loss')\n","    ax1.semilogy(epoch_xaxis,history.history['val_loss'],linewidth=2,label='Validation Loss')\n","    ax1.set_xlabel('Epoch [-]')\n","    ax1.set_ylabel('MSE [-]')\n","    ax1.grid()\n","    ax1.legend()\n","    ax2.plot(time, data_sets_norm[ch][u].iloc[SeedNo,1:],label ='Simulation',linewidth=2)\n","    ax2.plot(time,y_pred[0,:,:],label='Prediction',alpha=0.7,linewidth=2)\n","    ax2.set_ylabel('Scaled '+ch+' [-]')\n","    ax2.grid()\n","    ax2.legend()\n","    ax2.set_xlabel('Time [s]')\n","    a = ax3.hist(data_sets_norm[ch][u].iloc[SeedNo,1:], density=True, bins =20, orientation='horizontal',label ='Simulation',alpha = 0.75)\n","    _ = ax3.hist(y_pred[0,:,0],label='Prediction',bins = a[1], orientation='horizontal',alpha = 0.5,density=True)\n","    ax3.grid()\n","    ax3.legend()\n","    ax3.set_xlabel('Probablity [-]')\n","    fig.savefig(timestr+'_TCN_EncDec_UniWind_'+u+'_'+'Norm'+ch+'.png', dpi=400, facecolor='w',\n","             edgecolor='w', orientation='portrait', format='png', transparent=False)\n","  \n","\n","\n","\n","\n","\n"]},{"cell_type":"markdown","source":["# Testing: The hyper-parameter tuning/testing rig\n","\n","This piece of code is used to generate plots, while the hyper parameters in the `tcnae.py` code are changed. This was necessary to find an acceptable level of accuracy and keeping the track of what I was doing. The trainging part is not here, as it is a repeat of the loop above.\n","\n"],"metadata":{"id":"rZ1fjQgzUpjh"}},{"cell_type":"code","execution_count":null,"metadata":{"id":"17e66778-8437-4b04-81fa-4cf18b85acaf"},"outputs":[],"source":["timestr = tm.strftime('%Y%m%d%H%M%S')\n","time = np.linspace(0,600,12000)\n","seeds = np.arange(0,90)\n","rem_seeds = seeds[~np.in1d(np.arange(seeds.size),rnd_seed)]\n","SeedNo=random.choice(rem_seeds)\n","#rem_seeds = a[~np.in1d(np.arange(a.size),r)]\n","\n","#ch = 'MyTB'\n","#u = '10mps'\n","\n","\n","d = { 'a': str(tcn_ae.dilations),\n","      'b': str(tcn_ae.nb_filters),\n","      'c': str(tcn_ae.kernel_size),\n","      'd': str(tcn_ae.filters_conv1d),\n","      'e': str(tcn_ae.activation_conv1d),\n","      'f': str(tcn_ae.latent_sample_rate),\n","      'g': str(tcn_ae.tcn_activation),\n","      'h': str(tcn_ae.nb_stacks)}\n","\n","s = \"dilations={a}\\n nb_filters={b}\\n kernel_size={c}\\n filters_conv1d={d}\\n activation_conv1d={e}\\n latent_sample_rate={f}\\n tcn_activation={g}\\n nbs={h}\"\n","s.format(**d)\n","\n","\n","y_pred = tcn_ae.model.predict(data_sets_norm['wspdX'][u].iloc[SeedNo:SeedNo+1,1:])\n","#y_pred = loaded_model.predict(data_sets_norm['wspdX'][u].iloc[SeedNo:SeedNo+1,1:])\n","\n","print(y_pred.shape)\n","\n","fig,ax = plt.subplots(2,2,figsize=(20,14))\n","\n","ax=ax.flatten()\n","\n","ax[0].plot(time,data_sets_norm['wspdX'][u].iloc[SeedNo,1:])\n","ax[0].grid()\n","ax[0].set_ylabel('Normalized uniform wind speed [-]')\n","\n","ax[1].hist(data_sets_norm['wspdX'][u].iloc[SeedNo,1:],bins=30,density=True, orientation=u'horizontal')\n","ax[1].grid()\n","\n","\n","ax[2].plot(time, data_sets_norm[ch][u].iloc[SeedNo,1:],label ='Normalized Simulation')\n","ax[2].plot(time,y_pred[0,:,:],label='Normalized Prediction',alpha=0.6)\n","ax[2].set_ylabel('Normalized My @ Tower Bottom [-]')\n","ax[2].grid()\n","ax[2].legend()\n","ax[2].set_xlabel('Time [s]')\n","\n","a = ax[3].hist(data_sets_norm[ch][u].iloc[SeedNo,1:], density=True, bins =20, orientation='horizontal',label ='Normalized Simulation',alpha = 0.75)\n","_ = ax[3].hist(y_pred[0,:,0],label='Prediction',bins = a[1], orientation='horizontal',alpha = 0.5,density=True)\n","ax[3].grid()\n","ax[3].legend()\n","ax[3].set_xlabel('Probablity [-]')\n","plt.text(0.5,0.5,s.format(**d))\n","fig.savefig(timestr+'_TCN_EncDec_UniWind_'+u+'_'+'Norm'+ch+'.png', dpi=400, facecolor='w',\n","             edgecolor='w', orientation='portrait', format='png', transparent=False)"]}]}