Predictive Maintenance: [2-1] C-Mapss 데이터 분석편 (주제: Fault Classification)

심각하게 게으론 자의 C-Mapss 분석 2편이다.

금번 자료의 주제는 Fault Classification(혹은 Fault Diagnosis)이다.

 

C-Mapss 데이터의 목적은 RUL(Remaining Useful Lifetime) 예측이다. 하지만 근래 학계의 연구 동향을 봐도 RUL에 대해 진지하게 연구하고 노력하는 것들은 보기 어렵다(배터리 RUL제외). 왜냐? 어렵다. 미래 일을 어떻게 맞춰... 2022년 참석했던 PHM 학회에서도 현시대의 기술로는 RUL은 단순히 reference정도로 봐야한다는 평이 있었다. 나중에, Digital Twin이 완벽하게 구현된 시대에서는 RUL이 의미가 있을수는 있겠지만..

 

여하튼, 나도 위 의견에 매우-적극-굉장히 동의하므로 RUL 관점으로 분석은 하지 않을 것이다.

 

현실 세계에서는 설비의 고장 진단 즉, 고장 유형을 인식하는 것이 굉장히 중요하다.

 

아래 그림을 보자 (예지보전을 한다면 누구나 아는 그림)

제조 공정 중에 설비 고장은 해결해야하는 이슈다. 설비 자체에 대한 손실 뿐만아니라 제조 도메인에 따라서는 제조 고정 위에있는 모든 재료들을 버려야하기 때문이다. 특정 도메인의 경우에는 단위가 몇백억까지 가는 경우도 있다.

 

사람도 나이가 들면서 신체에 이상신호가 나타는 것이 당연하듯, 설비도 노후화 되면서 고장이 발생하는 것은 자연적인 현상이다

그렇기에 설비의 고장(=이상신호)을 조기에 탐지하고 설비 엔지니어가 정확한 정비 조치를 할 수 있도록 도와주는 것이 필요하다.

 

나는 "설비 엔지니어가 정확한 정비 조치"이 부분에 집중했다. 엔지니어에게 정확한 정비 조치를 할 수 있도록 도움을 줄 수 있는 것은 무엇일까 고민했다. 아직 자비스는 구현하지 못하는 실력이라, 고장 유형을 정확하게 알려주기만해도 정비 조치에는 큰 도움이 될 수 있을 것으로 생각했다.

 

 

Abstract가 너무 길었다. 요약하자면,

 

- 설비의 이상 상태를 정확하게 파악하여 엔지니어의 정비 조치를 돕겠다

 

- 필자는 이상 상태를 정확하게 파악하는 것을 "고장 유형 분류"라고 정의하고, 연구를 하겠다.

---

 

1. Fault Classification 관점의 EDA

 

금번 포스팅 자료의 target dataset은 C-Mapss의 FD001 train 데이터이다.

 

data description을 뒤져보면, FD001의 고장 유형은 한가지라고 나온다(1편 참조).

 

진짜 그럴까?

 

먼저 데이터를 파헤쳐보자.각종 자잘한 전처리는 1편을 참조하시라. 이전 편을 통해 분석 대상 파라미터는 정리가 끝났다.

- 코드는 아래 참고

path = '../1.data/cmapss/train_FD001.txt'

# define column names for easy indexing
index_names = ['key', 'time_cycles']
setting_names = ['setting_1', 'setting_2', 'setting_3']
sensor_names = ['s_{}'.format(i) for i in range(1,22)] 
col_names = index_names + setting_names + sensor_names

# read data
sf = pd.read_csv((path), sep='\s+', header=None, names=col_names)

delete_columns = ['time_cycles', 'setting_1', 'setting_2', 'setting_3', 's_1', 's_5', 's_6', 's_10', 's_16', 's_18', 's_19']
for s in delete_columns:
    del sf[s]
sf

 

FD001은 unit_nr라는 파라미터를 가지고 있다. 각 엔진을 의미하는 것이고, 엔진이 가동된 시점부터 고장날때까지(called "Run to Failure)" 수집된 데이터다. FD001의 unit_nr은 1~100을 가지고 있으므로 100번의 run to failure 이력이 있는 것이다.

 

그렇다면, data description에 따라 100가지의 이력들이 모두 하나의 고장 유형을 가지고 있어야한다.

진짜 그런지 살펴보겠다.

 

1.1 Dimension Reduction (Visualization)

 

모든 이력들의 초반 부분은 엔진의 상태가 정상이라고 조건을 추가하겠다. 대략 초반 50 point정도만 정상이라고 간주하겠다.

  

추가로 코드에서 "key"는 실제 데이터의 "unit_nr"과 동일.

sf_start_end = pd.DataFrame()
for k in sf['key'].unique():
    df_key = sf[sf['key']==k]
    
    for col in df_key.columns.tolist()[1:]:
        val = df_key[col].values
        val = (val-val[:50].mean())/val[:50].std()
        df_key[col] = val
        
    idx      = df_key.index.tolist()
    all_data = df_key
    sf_start_end = pd.concat([sf_start_end, all_data])
    
    
sf_start_end = sf_start_end.reset_index(drop=True)

keys = sf_start_end['key']
del sf_start_end['key']

 

이중 for문의 경우 데이터 전처리를 위한 것이고, 그 아래에 있는 all_data, sf_start_end는 Dimension reduction을 통한 시각화를 하기위해 만들었다.

 

이제 U-map을 통해 dimension reduction을 해보겠다.

import umap

reducer = umap.UMAP()
embedding = reducer.fit_transform(sf_start_end)

plt.figure(figsize=(16, 9))
plt.title('C-Mapss umap result')
for k in keys.unique():
    idx = keys[keys==k].index.tolist()
    plt.scatter(embedding[idx, 0], embedding[idx, 1], c=idx, cmap='jet', s=10, alpha=.3)

    trend_x = pd.Series(embedding[idx, 0]).rolling(window=10, min_periods=1).mean()
    trend_y = pd.Series(embedding[idx, 1]).rolling(window=10, min_periods=1).mean()
    plt.plot(trend_x, trend_y, c='k', alpha=.1)
plt.show()

 

umap을 사용한 이유는 t-sne는 실험을 할때마다 모델을 새로 생성해야한다는 단점이 있고, pca는 고차원에서는 umap에 비해 비교적 약한 성능을 보이는 것 같았다. 따라서 umap을 사용했다. 논리적인 근거는 아니지만, 본 포스팅에서 중요한 내용이 아니므로 건너간다.

 

umap을 돌린 결과, 정상에서 고장으로 갈때 3가지 방향으로 가는 것으로 보여진다.

 

이상하다. data description에서는 고장 유형이 한가지라고 했는데, 왜 3개의 방향으로 나뉘어질까? 하나의 방향으로 가야하는거 아닐까?

내가 이상한건지, description이 이상한건지 제대로 확인해보자.

 

정말 3가지의 방향으로 가는 것이 맞는지, 각 unit_nr의 끝에 있는 데이터 포인트 20개씩만 plotting 해본다.

embedding_result = pd.DataFrame()
plt.figure(figsize=(16, 9))
plt.title('C-Mapss umap result - end 20point')
for k in keys.unique():
    idx = keys[keys==k].index.tolist()
    # plt.scatter(embedding[idx, 0], embedding[idx, 1], c=idx, cmap='jet', s=10, alpha=.3)
    plt.scatter(embedding[idx[-20:], 0], embedding[idx[-20:], 1], c='C0', s=10, alpha=.6)
    
    embedding_result = pd.concat([embedding_result, pd.DataFrame(embedding[idx[-20:], :])])

plt.show()

위 그림을 보면 데이터의 끝 부분들이 3개 정도의 영역을 가지고 군집되어 있는 것으로 보여진다.

 

 

1.2 Clustering by K-means

각 unit_nr들의 고장 유형이 하나라면, 각 군집별로 파라미터를 plotting 했을 때 패턴이 유사할 것이다.

따라서 위 데이터를 기준으로 잡고 각 군집별로 label을 부여해서 plotting 결과를 보자.

label 부여 방법은 k-means를 사용(k=3)할 거다.

 

from sklearn.cluster import KMeans
embedding_result = embedding_result.reset_index(drop=True

kmeans = KMeans(n_clusters=3, n_init="auto").fit(embedding_result)

embedding_result['keys'] = 0
for i, st in enumerate(range(0, 2000, 20)):
    en = st + 20
    embedding_result['keys'].loc[st:en-1] = i + 1
    
embedding_result['label'] = kmeans.labels_

plt.figure(figsize=(16, 9))
plt.title('C-Mapss umap result')

for label in np.unique(kmeans.labels_):
    idx_label = np.where(kmeans.labels_==label)[0]
    plt.scatter(embedding_result.loc[idx_label, 0], embedding_result.loc[idx_label, 1], c='C%d'%(int(label)), s=10, alpha=.5, label=label)
    
plt.legend()
plt.show()

원했던 대로 3개의 군집으로 잘 나뉘어 졌다.

 

** 2024 07 08  내용 추가 **

하지만 K를 2로 설정하는 것이 정말 맞는 방법일까? K-means의 결과에 대해 평가할 수 있는 Silhouette index 기법을 사용한다.

Silhouette index란, 클러스터링(군집화) 결과를 평가하는 기법으로서, 각 데이터마다 속한 군집과 유사한 군집간의 거리를 계산하여 평가하는 기법이다. 해당 알고리즘에 대한 자세한 포스팅은 추후(아마도 안함...) 진행하겠다.

from sklearn.metrics import silhouette_score

result = []
for k in range(2, 15):
   X = embedding_result.iloc[:, :2]
   kmeans = KMeans(n_clusters=k, n_init="auto").fit(X)
   labels = kmeans.labels_
   result.append([k, silhouette_score(X, labels, metric='euclidean')])
result = np.array(result)

plt.figure(figsize=(16, 6))
plt.title('Silhouette Score Result about C-MAPSS')
plt.plot(result[:, 0], result[:, 1], marker='.', alpha=.6)
plt.show()

Silhouette index 결과를 봤을 때, K가 2일 경우에 silhouette index의 값이 가장 높게 나타났다.

따라서 K가 2일 경우에 가장 군집의 결과가 좋다고 평가할 수 있으며, 이후 분석에서 K는 2로 설정하여 진행하겠다.

 

 

1.3 Clustering Result Analysis

이제 전체 데이터의 parameter를 plotting해서 각 군집들 간의 차이가 얼마나 있는지 살펴보겠다.

key_cluster_reuslt = {}
for k in embedding_result['keys'].unique():
    label_result = embedding_result[embedding_result['keys']==k]['label']
    # print(label_result.value_counts().index.tolist()[0], label_result.value_counts())
    key_cluster_reuslt[k] = label_result.value_counts().index.tolist()[0]
    
    
import seaborn as sns

fig, ax = plt.subplots(5, 3, figsize=(40, 25))
for key in sf['key'].unique():
    color = 'C%d'%(key_cluster_reuslt[key])
    
    df_key = sf[sf['key']==key]
    for i, col in enumerate(df_key.columns.tolist()[1:]):
        row_ = i % 5
        col_ = i // 5
        ax[row_, col_].set_title(col)

        idx = np.arange(0, len(df_key))
        sns.lineplot(   ax=ax[row_, col_], x=idx, y=df_key[col], color=color, alpha=.3)
        sns.scatterplot(ax=ax[row_, col_], x=idx, y=df_key[col], color=color)

ax[0, 1].axvspan(0, 365, color='r', alpha=.1)
ax[4, 1].axvspan(0, 365, color='r', alpha=.1)

plt.show()
plt.close()

각 군집들에 대해 plotting한 결과다. 대체로 각 군집간의 차이는 없는 것으로 보인다. 각 엔진의 시작부터 끝까지 패턴은 거의 유사하다.

하지만 s_9, s_14 파라미터는 군집간의 패턴 차이가 보인다. 시작은 유사하나 끝에서 패턴이 달라지는 것으로 보여진다. 

 

중요한 parameter(s_9, s_14)만 plotting해보겠다.

fig, ax = plt.subplots(1, 2, figsize=(20, 6))

for key in sf['key'].unique():
    color = 'C%d'%(key_cluster_reuslt[key])

    df_key = sf[sf['key']==key]
    for i, col in enumerate(['s_9', 's_14']):
        if col == 's_9':
            row_ = 0
        else:
            row_ = 1
        ax[row_].set_title(col)

        idx = np.arange(0, len(df_key))
        sns.lineplot(   ax=ax[row_], x=idx, y=df_key[col], color=color, alpha=.3)
        sns.scatterplot(ax=ax[row_], x=idx, y=df_key[col], color=color)

ax[0].axvspan(0, 365, color='r', alpha=.1)
ax[1].axvspan(0, 365, color='r', alpha=.1)

plt.show()
plt.close()

대부분의 파라미터들은 군집별로 큰 차이가 보이지 않았다. 하지만 s_9, 14 파라미터의 경우 각 군집별로 패턴이 다르다.

결과적으로, 엔진이 고장날 때 패턴이 2가지 유형의 변화를 나타내며 끝났다는 것이다. 따라서 FD001의 고장 유형은 2가지라고 표현할 수 있다.

 

 

 

이상으로 C-Mapss 데이터를 이용한 고장 유형 분류 1단계(like EDA? for fault classification)를 진행했다.

매우 간단하지만, 정확한 fault classification을 하기 위한 첫단계로 이해해주시면 되겠다.

 

1번 포스팅(현재)에서는 fault classification을 위한 EDA를 진행했다.

2번 포스팅에서는 3가지 고장 유형들에 라벨을 부여하고, 딥러닝 기반의 모델에 학습시켜 잘 분류하는지 살펴보겠다.

3번 포스팅에서는 실제 현실에서 생길 수 있는 조건(open set recognition)을 추가해서 문제를 잘 해결하는지 살펴보겠다.