RUL_prediction_using_1D_CNN.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "RUL prediction using 1D CNN.ipynb",
      "private_outputs": true,
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uvL9GgtViveO"
      },
      "outputs": [],
      "source": [
        ""
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import tensorflow as tf\n",
        "from tensorflow.keras import layers\n",
        "from tensorflow.keras.models import Sequential\n",
        "\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import sklearn\n",
        "from sklearn.model_selection import train_test_split\n",
        "\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "from sklearn.metrics import mean_squared_error\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "np.random.seed(34)\n",
        "print(\"Tensorflow Version: \", tf.__version__)\n",
        "print(\"Numpy version: \", np.__version__)\n",
        "print(\"Pandas version: \", pd.__version__)\n",
        "print(\"Scikit-learn version: \", sklearn.__version__)"
      ],
      "metadata": {
        "id": "-FyRz9HCiwvl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "columns=[\"id\",\"cycle\",\"op1\",\"op2\",\"op3\",\"sensor1\",\"sensor2\",\"sensor3\",\"sensor4\",\"sensor5\",\"sensor6\",\"sensor7\",\"sensor8\",\n",
        "         \"sensor9\",\"sensor10\",\"sensor11\",\"sensor12\",\"sensor13\",\"sensor14\",\"sensor15\",\"sensor16\",\"sensor17\",\"sensor18\",\"sensor19\"\n",
        "         ,\"sensor20\",\"sensor21\"]"
      ],
      "metadata": {
        "id": "IGL0QKJDjATn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "train_data = pd.read_csv(\"/content/drive/MyDrive/train_FD001 (1).txt\", sep= \"\\s+\", header = None,names=columns )\n"
      ],
      "metadata": {
        "id": "P4pvwHLMjs8H"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "train_data. head() \n"
      ],
      "metadata": {
        "id": "qdvTgb-ZkZbr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def add_rul(g):\n",
        "    g['RUL'] = max(g['cycle']) - g['cycle']\n",
        "    return g\n",
        "\n",
        "train = train_data.groupby('id').apply(add_rul)\n",
        "train.head() "
      ],
      "metadata": {
        "id": "zPdkAA--kd2b"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "cnt_train=train_data[[\"id\", \"cycle\"]].groupby(\"id\").max().sort_values(by=\"cycle\", ascending=True)\n",
        "cnt_ind = [str(i) for i in cnt_train.index.to_list()]\n",
        "cnt_val = list(cnt_train.cycle.values)\n",
        "\n",
        "plt.style.use(\"fivethirtyeight\")\n",
        "fig, ax = plt.subplots(figsize=(12, 30))\n",
        "rects = ax.barh(cnt_ind, cnt_val)\n",
        "ax.set_title(\"Number of Cycles of Different Engines\", fontsize=22)\n",
        "ax.set_ylabel(\"Engine Id\", fontsize=18)\n",
        "ax.set_xlabel(\"Number of Cycles\", fontsize=18)\n",
        "ax.grid(True)\n",
        "\n",
        "for rect in rects:\n",
        "    width = rect.get_width()\n",
        "    plt.text(1.05*rect.get_width(), rect.get_y()+0.5*rect.get_height(),\n",
        "             '%d' % int(width),\n",
        "             ha='center', va='center')"
      ],
      "metadata": {
        "id": "P88WErtHksfQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "subset_stats = train_data.agg(['mean', 'std']).T[2:]\n",
        "ax = subset_stats.plot.bar(figsize=(18,10), \n",
        "                           y=\"std\")\n",
        "ax.set_title(\"Standard Deviation of Every Sensor From Every Engine\", fontsize=22)\n",
        "ax.set_xlabel(\"Sensor\", fontsize=18)\n",
        "ax.set_ylabel(\"Standard Deviation\", fontsize=18)\n",
        "\n",
        "for p in ax.patches:\n",
        "    ax.annotate(str(round(p.get_height(),3)), (p.get_x() * 1.005, p.get_height() * 1.005))\n"
      ],
      "metadata": {
        "id": "-GneI5qFkxog"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def process_targets(data_length, early_rul = None):\n",
        "    \"\"\" \n",
        "    Takes datalength and earlyrul as input and \n",
        "    creates target rul.\n",
        "    \"\"\"\n",
        "    if early_rul == None:\n",
        "        return np.arange(data_length-1, -1, -1)\n",
        "    else:\n",
        "        early_rul_duration = data_length - early_rul\n",
        "        if early_rul_duration <= 0:\n",
        "            return np.arange(data_length-1, -1, -1)\n",
        "        else:\n",
        "            return np.append(early_rul*np.ones(shape = (early_rul_duration,)), np.arange(early_rul-1, -1, -1))\n",
        "def process_input_data_with_targets(input_data, target_data = None, window_length = 1, shift = 1):\n",
        "    \"\"\"Depending on values of window_length and shift, this function generates batchs of data and targets \n",
        "    from input_data and target_data.\n",
        "    \n",
        "    Number of batches = np.floor((len(input_data) - window_length)/shift) + 1\n",
        "    \n",
        "    **We don't check input dimensions uisng exception handling. So readers should be careful while using these\n",
        "    functions. If input data are not of desired dimension, either error occurs or something undesirable is \n",
        "    produced as output.**\n",
        "    \n",
        "    Arguments:\n",
        "        input_data: input data to function (Must be 2 dimensional)\n",
        "        target_data: input rul values (Must be 1D array)s\n",
        "        window_length: window length of data\n",
        "        shift: Distance by which the window moves for next batch. This is closely related to overlap\n",
        "               between data. For example, if window length is 30 and shift is 1, there is an overlap of \n",
        "               29 data points between two consecutive batches.\n",
        "        \n",
        "    \"\"\"\n",
        "    num_batches = np.int(np.floor((len(input_data) - window_length)/shift)) + 1\n",
        "    num_features = input_data.shape[1]\n",
        "    output_data = np.repeat(np.nan, repeats = num_batches * window_length * num_features).reshape(num_batches, window_length,\n",
        "                                                                                                  num_features)\n",
        "    if target_data is None:\n",
        "        for batch in range(num_batches):\n",
        "            output_data[batch,:,:] = input_data[(0+shift*batch):(0+shift*batch+window_length),:]\n",
        "        return output_data\n",
        "    else:\n",
        "        output_targets = np.repeat(np.nan, repeats = num_batches)\n",
        "        for batch in range(num_batches):\n",
        "            output_data[batch,:,:] = input_data[(0+shift*batch):(0+shift*batch+window_length),:]\n",
        "            output_targets[batch] = target_data[(shift*batch + (window_length-1))]\n",
        "        return output_data, output_targets\n",
        "def process_test_data(test_data_for_an_engine, window_length, shift, num_test_windows = 1):\n",
        "    \"\"\" This function takes test data for an engine as first input. The next two inputs\n",
        "    window_length and shift are same as other functins. \n",
        "    \n",
        "    Finally it takes num_test_windows as the last input. num_test_windows sets how many examplles we\n",
        "    want from test data (from last). By default it extracts only the last example.\n",
        "    \n",
        "    The function return last examples and number of last examples (a scaler) as output. \n",
        "    We need the second output later. If we are extracting more than 1 last examples, we have to \n",
        "    average their prediction results. The second scaler halps us do just that.\n",
        "    \"\"\"\n",
        "    max_num_test_batches = np.int(np.floor((len(test_data_for_an_engine) - window_length)/shift)) + 1\n",
        "    if max_num_test_batches < num_test_windows:\n",
        "        required_len = (max_num_test_batches -1)* shift + window_length\n",
        "        batched_test_data_for_an_engine = process_input_data_with_targets(test_data_for_an_engine[-required_len:, :],\n",
        "                                                                          target_data = None,\n",
        "                                                                          window_length = window_length, shift = shift)\n",
        "        return batched_test_data_for_an_engine, max_num_test_batches\n",
        "    else:\n",
        "        required_len = (num_test_windows - 1) * shift + window_length\n",
        "        batched_test_data_for_an_engine = process_input_data_with_targets(test_data_for_an_engine[-required_len:, :],\n",
        "                                                                          target_data = None,\n",
        "                                                                          window_length = window_length, shift = shift)\n",
        "        return batched_test_data_for_an_engine, num_test_windows"
      ],
      "metadata": {
        "id": "bf3ug2kKlAmD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "bIOW9k7xl4C4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "test_data = pd.read_csv(\"/content/drive/MyDrive/test_FD001.txt\", sep = \"\\s+\", header = None,names=columns)\n",
        "true_rul = pd.read_csv(\"/content/drive/MyDrive/RUL_FD001.txt\", sep = '\\s+', header = None)\n",
        "\n",
        "window_length = 30\n",
        "shift = 1\n",
        "early_rul = 125            \n",
        "processed_train_data = []\n",
        "processed_train_targets = []\n",
        "\n",
        "# How many test windows to take for each engine. If set to 1 (this is the default), only last window of test data for \n",
        "# each engine is taken. If set to a different number, that many windows from last are taken. \n",
        "# Final output is the average output of all windows.\n",
        "num_test_windows = 5     \n",
        "processed_test_data = []\n",
        "num_test_windows_list = []\n",
        "\n",
        "#columns_to_be_dropped = [0,1,2,3,4,5,9,10,14,20,22,23]\n",
        "columns_to_be_dropped =['id','op1','op2','op3', 'sensor1', 'sensor5','sensor6','sensor7', 'sensor10',\n",
        "                 'sensor16', 'sensor18', 'sensor19']\n",
        "train_data_first_column = train_data['id']\n",
        "test_data_first_column = test_data['id']\n",
        "\n",
        "# Scale data for all engines\n",
        "scaler = MinMaxScaler(feature_range = (-1,1))\n",
        "train_data = scaler.fit_transform(train_data.drop(columns = columns_to_be_dropped))\n",
        "test_data = scaler.transform(test_data.drop(columns = columns_to_be_dropped))\n",
        "\n",
        "train_data = pd.DataFrame(data = np.c_[train_data_first_column, train_data])\n",
        "test_data = pd.DataFrame(data = np.c_[test_data_first_column, test_data])\n",
        "\n",
        "num_train_machines = len(train_data[0].unique())\n",
        "num_test_machines = len(test_data[0].unique())\n",
        "\n",
        "# Process training and test data sepeartely as number of engines in training and test set may be different.\n",
        "# As we are doing scaling for full dataset, we are not bothered by different number of engines in training and test set.\n",
        "\n",
        "# Process trianing data\n",
        "for i in np.arange(1, num_train_machines + 1):\n",
        "    temp_train_data = train_data[train_data[0] == i].drop(columns = [0]).values\n",
        "    \n",
        "    # Verify if data of given window length can be extracted from training data\n",
        "    if (len(temp_train_data) < window_length):\n",
        "        print(\"Train engine {} doesn't have enough data for window_length of {}\".format(i, window_length))\n",
        "        raise AssertionError(\"Window length is larger than number of data points for some engines. \"\n",
        "                             \"Try decreasing window length.\")\n",
        "        \n",
        "    temp_train_targets = process_targets(data_length = temp_train_data.shape[0], early_rul = early_rul)\n",
        "    data_for_a_machine, targets_for_a_machine = process_input_data_with_targets(temp_train_data, temp_train_targets, \n",
        "                                                                                window_length = window_length, shift = shift)\n",
        "    \n",
        "    processed_train_data.append(data_for_a_machine)\n",
        "    processed_train_targets.append(targets_for_a_machine)\n",
        "\n",
        "processed_train_data = np.concatenate(processed_train_data)\n",
        "processed_train_targets = np.concatenate(processed_train_targets)\n",
        "\n",
        "# Process test data\n",
        "for i in np.arange(1, num_test_machines + 1):\n",
        "    temp_test_data = test_data[test_data[0] == i].drop(columns = [0]).values\n",
        "    \n",
        "    # Verify if data of given window length can be extracted from test data\n",
        "    if (len(temp_test_data) < window_length):\n",
        "        print(\"Test engine {} doesn't have enough data for window_length of {}\".format(i, window_length))\n",
        "        raise AssertionError(\"Window length is larger than number of data points for some engines. \"\n",
        "                             \"Try decreasing window length.\")\n",
        "    \n",
        "    # Prepare test data\n",
        "    test_data_for_an_engine, num_windows = process_test_data(temp_test_data, window_length = window_length, shift = shift,\n",
        "                                                             num_test_windows = num_test_windows)\n",
        "    \n",
        "    processed_test_data.append(test_data_for_an_engine)\n",
        "    num_test_windows_list.append(num_windows)\n",
        "\n",
        "processed_test_data = np.concatenate(processed_test_data)\n",
        "true_rul = true_rul[0].values\n",
        "\n",
        "# Shuffle training data\n",
        "index = np.random.permutation(len(processed_train_targets))\n",
        "processed_train_data, processed_train_targets = processed_train_data[index], processed_train_targets[index]\n",
        "\n",
        "print(\"Processed trianing data shape: \", processed_train_data.shape)\n",
        "print(\"Processed training ruls shape: \", processed_train_targets.shape)\n",
        "print(\"Processed test data shape: \", processed_test_data.shape)\n",
        "print(\"True RUL shape: \", true_rul.shape)"
      ],
      "metadata": {
        "id": "0C5T70HelcCn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "processed_train_data, processed_val_data, processed_train_targets, processed_val_targets = train_test_split(processed_train_data,\n",
        "                                                                                                            processed_train_targets,\n",
        "                                                                                                            test_size = 0.2,\n",
        "                                                                                                            random_state = 83)\n",
        "print(\"Processed train data shape: \", processed_train_data.shape)\n",
        "print(\"Processed validation data shape: \", processed_val_data.shape)\n",
        "print(\"Processed train targets shape: \", processed_train_targets.shape)\n",
        "print(\"Processed validation targets shape: \", processed_val_targets.shape)"
      ],
      "metadata": {
        "id": "GxwuqdTBnOpX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def create_compiled_model():\n",
        "    model = Sequential([\n",
        "        layers.Conv1D(256, 7, activation = \"relu\", input_shape = (window_length, processed_train_data.shape[2])),\n",
        "        layers.Conv1D(96, 7, activation = \"relu\"),\n",
        "        layers.Conv1D(32, 7, activation = \"relu\"),\n",
        "        layers.GlobalAveragePooling1D(),\n",
        "        layers.Dense(64, activation = \"relu\"),\n",
        "        layers.Dense(128, activation = \"relu\"),\n",
        "        layers.Dense(1)\n",
        "    ])\n",
        "    model.compile(loss = \"mse\", optimizer = tf.keras.optimizers.Adam(learning_rate=0.001))\n",
        "    return model"
      ],
      "metadata": {
        "id": "nC-V6_8onS-m"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def scheduler(epoch):\n",
        "    if epoch < 10:\n",
        "        return 0.001\n",
        "    else:\n",
        "        return 0.0001\n",
        "callback = tf.keras.callbacks.LearningRateScheduler(scheduler, verbose = 1)\n",
        "model = create_compiled_model()\n",
        "history = model.fit(processed_train_data, processed_train_targets, epochs = 30,\n",
        "                    validation_data = (processed_val_data, processed_val_targets),\n",
        "                    callbacks = callback,\n",
        "                    batch_size = 64, verbose = 2)"
      ],
      "metadata": {
        "id": "YoFMUcEinbQY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "rul_pred = model.predict(processed_test_data).reshape(-1)\n",
        "preds_for_each_engine = np.split(rul_pred, np.cumsum(num_test_windows_list)[:-1])\n",
        "mean_pred_for_each_engine = [np.average(ruls_for_each_engine, weights = np.repeat(1/num_windows, num_windows)) \n",
        "                             for ruls_for_each_engine, num_windows in zip(preds_for_each_engine, num_test_windows_list)]\n",
        "RMSE = np.sqrt(mean_squared_error(true_rul, mean_pred_for_each_engine))\n",
        "print(\"RMSE: \", RMSE)"
      ],
      "metadata": {
        "id": "dJhYkb-5nrjD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "tf.keras.models.save_model(model, \"FD001_1D_CNN_piecewise_RMSE_\"+ str(np.round(RMSE, 4)) + \".h5\")"
      ],
      "metadata": {
        "id": "tuxSPI2kn0Vr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "indices_of_last_examples = np.cumsum(num_test_windows_list) - 1\n",
        "preds_for_last_example = np.concatenate(preds_for_each_engine)[indices_of_last_examples]\n",
        "\n",
        "RMSE_new = np.sqrt(mean_squared_error(true_rul, preds_for_last_example))\n",
        "print(\"RMSE (Taking only last examples): \", RMSE_new)"
      ],
      "metadata": {
        "id": "4uPXQMGCn7Jt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def compute_s_score(rul_true, rul_pred):\n",
        "    \"\"\"\n",
        "    Both rul_true and rul_pred should be 1D numpy arrays.\n",
        "    \"\"\"\n",
        "    diff = rul_pred - rul_true\n",
        "    return np.sum(np.where(diff < 0, np.exp(-diff/13)-1, np.exp(diff/10)-1))\n",
        "s_score = compute_s_score(true_rul, preds_for_last_example)\n",
        "print(\"S-score: \", s_score)"
      ],
      "metadata": {
        "id": "keWdL7HyoEtM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot true and predicted RUL values\n",
        "plt.plot(true_rul, label = \"True RUL\", color = \"red\")\n",
        "plt.plot(preds_for_last_example, label = \"Pred RUL\", color = \"blue\")\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "1Q2sB68loOjH"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}