NASA Turbofan Engine Datasets : 教師なし学習

# 初めに

PHM2008のNASA Turbofan Engineのデータセットを利用して分析してみます。
https://data.nasa.gov/Aerospace/CMAPSS-Jet-Engine-Simulated-Data/ff5v-kuh6
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/02229b9b-300d-1fdd-ea8d-974a93b4deaa.png)

# データセットの準備

## データセットの説明（英文)
The data are provided as a zip-compressed text file with 26 columns of numbers, separated by spaces. Each row is a snapshot of data taken during a single operational cycle, each column is a different variable.

The columns correspond to:

unit number
time, in cycles
operational setting 1
operational setting 2
operational setting 3
sensor measurement 1
sensor measurement 2 ...
sensor measurement 26
Data Set: FD001

Train trjectories: 100
Test trajectories: 100
Conditions: ONE (Sea Level)
Fault Modes: ONE (HPC Degradation)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/b9623cec-a879-5baf-43a2-c72280319710.png)

# 探索的データ解析

```python:
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
%matplotlib inline

df_train = pd.read_csv('/PHM-IEEE-2008/train_FD001.csv', header = None)
```

```python:
columns_train = ['unit_ID','cycles','setting_1','setting_2','setting_3','T2','T24','T30','T50','P2','P15','P30','Nf',
           'Nc','epr','Ps30','phi','NRf','NRc','BPR','farB','htBleed','Nf_dmd','PCNfR_dmd','W31','W32' ]

df_train.columns = columns_train
df_train
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/0eff6213-c66a-2ec1-1db6-c85fe4ad07ca.png)

## plot
```python:
fig, axes = plt.subplots(6,5, figsize=(30,30))

ax = axes.ravel()

for i in range(len(df_train.columns)):
    
    ax[i].scatter(df_train['cycles'], df_train.iloc[:,i], color='blue')
    ax[i].set_title(df_train.columns[i])
    ax[i].set_yticks(())
    ax[i].grid()

fig.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/7d3af808-3ed0-0c00-79f0-77d8a4e8b94b.png)


# データの前処理

## 移動平均(n=5)
ノイズ除去

```python:
#making new dataframe
df_train_sma = pd.DataFrame(columns=columns_train)

df_train.columns

columns_to_move_average = ['setting_1', 'setting_2', 'setting_3', 'T2', 'T24',
       'T30', 'T50', 'P2', 'P15', 'P30', 'Nf', 'Nc', 'epr', 'Ps30', 'phi',
       'NRf', 'NRc', 'BPR', 'farB', 'htBleed', 'Nf_dmd', 'PCNfR_dmd', 'W31',
       'W32']

for name in columns_to_move_average :
    df_train_sma[name] = df_train[name].rolling(5).mean()

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/3e711787-9fd1-c950-bdb4-77803ac27286.png)

unit_IDとcyclesのデータを戻す

```python:
df_train_sma['unit_ID'] = df_train['unit_ID']
df_train_sma['cycles'] = df_train['cycles']
df_train_sma
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/fb9beac5-a3d1-a6d4-a295-9844b87b9f55.png)

NAをdrop

```python:
# drop NA
df_train_sma.dropna(inplace=True)
df_train_sma

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/7b4417f7-a9f7-9375-2fca-069044781a29.png)

indexのreset
```python:
df_train_sma.reset_index(drop=True, inplace=True)
df_train_sma

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/715648c7-cc6e-0656-a901-ca4a472c3322.png)

再びplot
```python:
fig, axes = plt.subplots(6,5, figsize=(30,30))

ax = axes.ravel()

for i in range(len(df_train.columns)):
    
    ax[i].scatter(df_train_sma['cycles'], df_train_sma.iloc[:,i], color='red')
    ax[i].set_title(df_train.columns[i])
    ax[i].set_yticks(())
    ax[i].grid()

fig.tight_layout()
plt.show()

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/12bf5d63-523c-dfbb-77c1-26d8571dc16d.png)


## Feature選択

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/58bba582-88d7-3818-cd60-f2ee1e7fc8b6.png)

columns_to_drop_train = ['setting_1', 'setting_2', 'setting_3', 'T2','P2', 'P15', 'epr','farB', 'Nf_dmd', 'PCNfR_dmd', 'W31','W32']

```python:
df_train.columns
columns_to_drop_train = [ 'setting_1', 'setting_2', 'setting_3', 
                          'T2','P2', 'P15', 'epr','farB', 'Nf_dmd', 'PCNfR_dmd',
                          'W31','W32']
df_train_feature = df_train_sma.drop(columns=columns_to_drop_train)
df_train_feature

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/f7fa287c-a2bc-b023-0627-0934260080b8.png)




## Standarization

Scalingを行う。unit_ID、cyclesはscaling対象ではないことに気を付けましょう。
```python:
# Scaling
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
df_pca = df_train_feature.drop(columns=['unit_ID', 'cycles'])

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/2ec5a80d-8e34-958d-50df-c234e10f14a8.png)


```python:
scaler.fit(df_pca)
df_scaled = pd.DataFrame(scaler.transform(df_pca), columns = df_pca.columns)
df_scaled

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/79bccfc0-fe5b-ab17-538b-20b3a2726c01.png)

# データ分析（PCA）

```python:
#PCA
from sklearn.decomposition import PCA
pca = PCA(n_components=4)
pca.fit(df_scaled)
pca_result = pca.transform(df_scaled)

```
```python:
print(f'original data: {df_scaled.shape}')
print(f'after PCA data: {pca_result.shape}')
print(f' Explaiend variance ratio : {pca.explained_variance_ratio_}')

```
PCA1 77%とPCA2 18% 合計で95% 悪くないですね。
```
original data: (20627, 12)
after PCA data: (20627, 4)
 Explaiend variance ratio : [0.76605196 0.18037961 0.01050908 0.00907406]

```


# データ分析結果の解析

pca結果をdataframe化する。その後の可視化、解析のため。
```python:
# pca result to dataframe
df_pca_result = pd.DataFrame(pca_result, columns=['PC1', 'PC2','PC3','PC4',])
df_pca_result
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/e9169a92-dcc6-19f4-bb82-0e40ddddac47.png)

## とにかく可視化

```python:
fig , ax = plt.subplots(1,1, figsize = (10,10))

ax.scatter(x=df_pca_result['PC1'], y=df_pca_result['PC2'])
ax.set_title('PCA')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')

ax.grid()

fig.tight_layout()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/7ed39ad1-2461-0f1e-ffd3-836b9ee5c478.png)

## cyclesの情報を入れてみよう。

```python:
df_pca_result_plus_cycles = pd.concat([df_pca_result, df_train_feature[['unit_ID', 'cycles']] ], axis=1)
df_pca_result_plus_cycles
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/bc10deff-a8a4-768c-ea1c-c7d0076a814b.png)

```python:
fig , ax = plt.subplots(1,1, figsize = (12,10))

mappable=ax.scatter(x=df_pca_result_plus_cycles['PC1'], y=df_pca_result_plus_cycles['PC2'],
           c=df_pca_result_plus_cycles['cycles'])
ax.set_title('PCA')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')

ax.grid()

fig.tight_layout()
fig.colorbar(mappable)
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/e141f28b-9408-b9e2-0a2c-3e2afdc9770a.png)

## まだよく分からないので、四つのエンジンだけ取り出してみる。

```python:
# Unit Number 35, 57, 79, 91
unit_list = [35, 57, 79, 91]
#unit_list = [1, 2, 3, 100]
fig , ax = plt.subplots(1,1, figsize = (10,10))

for number in unit_list:


    ax.scatter(x=df_pca_result_plus_cycles[df_pca_result_plus_cycles['unit_ID']==number]['PC1'],
                          y=df_pca_result_plus_cycles[df_pca_result_plus_cycles['unit_ID']==number]['PC2'],
           )
    
ax.set_title('PCA')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')    
ax.grid()
ax.legend(unit_list)
fig.tight_layout()

plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/d089deef-ef9f-299c-b95b-57546c1d887a.png)

## 1つのエンジンだけにして、cyclesで色をつけてみる。
```python:
fig , ax = plt.subplots(1,1, figsize = (12,10))
number = 5

mappable=ax.scatter(x=df_pca_result_plus_cycles[df_pca_result_plus_cycles['unit_ID']==number]['PC1'],
                    y=df_pca_result_plus_cycles[df_pca_result_plus_cycles['unit_ID']==number]['PC2'],
                    c=df_pca_result_plus_cycles[df_pca_result_plus_cycles['unit_ID']==number]['cycles'])
ax.set_title('PCA')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')

ax.grid()
ax.legend([number,])
fig.tight_layout()
fig.colorbar(mappable, ax=None)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/a77c6d7a-6435-3b9f-d5fa-1877e2d105b8.png)

PCAでエンジンの劣化状態が分かりますね。


# 参考文献

* [CMAPSS Jet Engine Simulated Data](https://data.nasa.gov/Aerospace/CMAPSS-Jet-Engine-Simulated-Data/ff5v-kuh6)

* [https://www.kaggle.com/code/mochmunirulichsan/predictive-maintenance-scheduling-fd001](https://www.kaggle.com/code/mochmunirulichsan/predictive-maintenance-scheduling-fd001)

* [https://pub.towardsai.net/predicting-the-remaining-useful-life-of-turbofan-engine-f38a17391cac](https://pub.towardsai.net/predicting-the-remaining-useful-life-of-turbofan-engine-f38a17391cac)


* [Remaining useful life prognosis of turbofan engines based on deep feature extraction and fusion](https://www.nature.com/articles/s41598-022-10191-2?fbclid=IwAR2BKQte_6YKFJ0hGesOo9jF0GfokayI79bW_7ivru-rdBlcGLus0hcW1U8)

* [Prediction of Remaining Useful Life of an Engine Based on Sensors; Building A Random Forest in Python](https://medium.com/@polanitzer/prediction-of-remaining-useful-life-of-an-engine-based-on-sensors-building-a-random-forest-in-ffad82c8a1c6)

* [Predictive Maintenance (PdM) of Aircraft Engine](https://github.com/archd3sai/Predictive-Maintenance-of-Aircraft-Engine?fbclid=IwAR1LbsG4gPckuaCswDdDaLPjwKXPD7PLwJ-PXs5szwh_kanWMpa0B-qQasY)

*  [https://github.com/mathworks/JP_Predictive_Maintenance_MATLAB/blob/master/docs/UnsupervisedLive_JP.md](https://github.com/mathworks/JP_Predictive_Maintenance_MATLAB/blob/master/docs/UnsupervisedLive_JP.md)

* https://github.com/mathworks/pdm_demo_by_matlab

* Course2 Statement: NASA Turbofan PCA](https://www.kaggle.com/code/brjapon/course2-statement-nasa-turbofan-pca)

* https://github.com/mathworks/JP_Predictive_Maintenance_MATLAB/blob/master/docs/UnsupervisedLive_JP.md

* https://jp.mathworks.com/videos/predictive-maintenance-with-matlab-a-prognostics-case-study-121138.html