발전소에서 나온 데이터 분석해보기

'''
Acknowledgements

Source:

Pınar Tüfekci, Çorlu Faculty of Engineering, Namık Kemal University, TR-59860 Çorlu, Tekirdağ, Turkey
Email: ptufekci '@' nku.edu.tr

Heysem Kaya, Department of Computer Engineering, Boğaziçi University, TR-34342, Beşiktaş, İstanbul, Turkey
Email: heysem '@' boun.edu.tr

'''
print("")

Combined Cycle Power Plant Data Set

데이터분석, 시각화와 예측 모델 구축에 적용하기

이 데이터 셋은 발전소로부터 얻은 5년치 이상의 센서 데이터 입니다.

공장의 시간당 전기 에너지 출력 (EP)을 예측하기위한
시간별 평균 주변 온도 (AT),
주변 압력 (AP),
상대 습도 (RH) 및
배출 진공 (V)으로 구성됩니다.

복합 화력 발전소 (CCPP)는 가스 터빈 (GT), 증기 터빈 (ST) 및 열 회수 증기 발생기로 구성됩니다.

CCPP에서 전기는 한 주기로 합쳐진 가스 터빈 및 증기 터빈에 의해 생성되고 한 터빈에서 다른 터빈으로 전달됩니다.

Vacuum이 Steam Turbine으로부터 통합되어 영향을받는 동안, 다른 주변 변수 중 세 가지가 GT 성능에 영향을 미칩니다.

컬럼 설명

(데이터는 다음의 링크에서 다운받으실 수 있습니다)

• EP

  • Net hourly electrical energy output (EP) 420.26-495.76 MW -> 공장의 시간당 전기 에너지 출력

• AT

  • Ambient Temperature (AT) in the range 1.81°C and 37.11°C -> 시간별 평균 주변 온도

• AP

  • Ambient Pressure (AP) in the range 992.89-1033.30 milibar -> 시간별 평균 주변 압력

• RH

  • Relative Humidity (RH) in the range 25.56% to 100.16% -> 시간별 평균 주변 습도

• V

  • Exhaust Vacuum (V) in teh range 25.36-81.56 cm Hg -> 시간펼 평균 배출 진공도

• 참고사항

  • 변수들은 데이터 획득시 공장 주변에 위치한 다양한 센서로부터 도출된 값들의 평균으로 이루어져 있습니다.
  • 변수는 정규화없이 주어집니다.

# 파이썬의 데이터 분석 도구인 Pandas 를 불러옵니다.
# Pandas 를 쉽게 설명드리면, 파이썬으로 엑셀을 다룰 수 있는 도구라고 볼 수 있습니다.
# 이 도구를 앞으로는 pd라는 축약어로 사용하기로 합니다.
import pandas as pd

# matplotlib로 그래프를 그릴 때, 바로 화면에 보이게끔 만들어 줍니다.
%matplotlib inline

Load Dataset

모든 데이터 분석의 시작은 데이터를 읽어오는 것입니다.

파일의 경로를 지정하는 방법에 주의하셔야 합니다.

만일 read_csv를 실행할 때 (FileNotFoundError)라는 이름의 에러가 난다면 경로가 제대로 지정이 되지 않은 것입니다.

다음의 링크에서 경로를 지정하는 법을 다루고 있습니다.

탐색기에서 파일이 있는 경로까지 가서 경로주소를 복사한 다음, path 변수에 문자열로 넣어주세요.

원 표시 또는 역슬래시가 있다면 그것을 슬래시 ( / ) 로 바꿔주세요.
맨 마지막에도 슬래시를 붙여주세요.

예)

C:\Users\one\Desktop\Factory\blade\data 이라면,

C:/Users/one/Desktop/Factory/blade/data/ 로 바꿔서

path = 'C:/Users/one/Desktop/Factory/blade/data/' 이렇게 path 변수에 문자열로 넣어주세요

# 판다스의 read_excel로 Folds5x2_pp.xlsx 파일을 읽어옵니다.

# 이 데이터셋은 Sheet1 부터 Sheet5 까지 데이터가 나뉘어져 있습니다.
# 각 시트별로 따로따로 불러와 각기 다른 변수(Sheet1 ~ Sheet5)에 넣어줍니다.

path = 'C:/Users/1990y/Downloads/Factory/power plant/'

path = path + 'Folds5x2_pp.xlsx'
Sheet1 = pd.read_excel(path, sheet_name='Sheet1')
Sheet2 = pd.read_excel(path, sheet_name='Sheet2')
Sheet3 = pd.read_excel(path, sheet_name='Sheet3')
Sheet4 = pd.read_excel(path, sheet_name='Sheet4')
Sheet5 = pd.read_excel(path, sheet_name='Sheet5')

Combine Dataset

Sheet1 부터 Sheet5 까지를 하나의 데이터프레임으로 묶기

# year 컬럼을 새로 만들어, 1년차는 1, 2년차는 2, ... 5년차는 5를 넣어줍니다.
# 섞였을 때 서로 구분하기 위함입니다.

Sheet1['year'] = 1
Sheet2['year'] = 2
Sheet3['year'] = 3
Sheet4['year'] = 4
Sheet5['year'] = 5

# pd.concat() 안에 리스트 형태로, 묶고자 하는 데이터프레임들을 넣어주면 됩니다.
# 왼쪽(Sheet) 의 아래에 오른쪽(Sheet2, Sheet3, Sheet4, Sheet5) 이 차례대로 붙는다고 생각하시면 됩니다.
# 이어붙인 값을 = 을 통해 combined 라는 변수에 할당합니다.

combined = pd.concat([Sheet1, Sheet2, Sheet3, Sheet4, Sheet5])

combined.describe()

	AT	V	AP	RH	PE	year
count	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000
mean	19.651231	54.305804	1013.259078	73.308978	454.365009	3.000000
std	7.452162	12.707362	5.938535	14.599658	17.066281	1.414228
min	1.810000	25.360000	992.890000	25.560000	420.260000	1.000000
25%	13.510000	41.740000	1009.100000	63.327500	439.750000	2.000000
50%	20.345000	52.080000	1012.940000	74.975000	451.550000	3.000000
75%	25.720000	66.540000	1017.260000	84.830000	468.430000	4.000000
max	37.110000	81.560000	1033.300000	100.160000	495.760000	5.000000

# 기존에 있던 컬럼들의 상호 비교를 위해 모든 변수들과의 관계에 대한 산점도 행렬을 구합니다.
# from pandas.plotting import * 는 판다스 안의 산점도 행렬 함수를 불러오기 위함입니다.
# 산점도 행렬 함수 안에 들어가는 내용은 다음과 같습니다.

# combined -> 보려고 하는 모든 데이터가 있는 데이터프레임
# c -> color 의 약자로, 산점도 안에서 구분하고자 하는 정보, 지금은 combined['label'] 이고, 이는 새 것과 헌 것을 구분합니다.
# figsize -> 보려고 하는 이미지의 크기 정보

from pandas.plotting import *

features = ['AT', 'V', 'AP', 'RH']
scatter_matrix(combined[features], 
                c = combined['PE'], 
                figsize=(25, 25))

# 산점도 행렬은 이후 여러번 사용할 것이기 때문에, 반복을 줄이기 위해 함수로 만들어줍니다.

def scatter(data, label):
    scatter_matrix(data, c = label, figsize=(25, 25))

파생변수 생성

# 데이터 분석자의 상식과 배경지식을 동원하여, 기존 특징들로부터 도출할 수 있는 새로운 변수입니다.
# 현재 데이터셋의 변수들은 1시간 평균치의 물리량이므로, 
# 1시간 전으로부터 환경변수의 변화율을 구할 수 있습니다.

# .pct_change(1) 을 하여, 이전 row 와의 퍼센트 차이를 구할 수 있습니다.

target_dset = combined

target_dset['1d_AT'] = target_dset['AT'].pct_change(1)
target_dset['1d_V'] = target_dset['V'].pct_change(1)
target_dset['1d_AP'] = target_dset['AP'].pct_change(1)
target_dset['1d_RH'] = target_dset['RH'].pct_change(1)

# .fillna(0) 을 하는 이유는, 첫 번째 행(row) 의 경우 이전 값이 존재하지 않기 때문에 
# NaN 이 있으므로, 이를 채워주기 위함입니다.
combined_change = target_dset.fillna(0)

# 잘 연산되었는지 확인하기 위해 .head() 로 첫 5 번째 행들을 알아봅니다.
combined_change.head()

	AT	V	AP	RH	PE	year	1d_AT	1d_V	1d_AP	1d_RH
0	14.96	41.76	1024.07	73.17	463.26	1	0.000000	0.000000	0.000000	0.000000
1	25.18	62.96	1020.04	59.08	444.37	1	0.683155	0.507663	-0.003935	-0.192565
2	5.11	39.40	1012.16	92.14	488.56	1	-0.797061	-0.374206	-0.007725	0.559580
3	20.86	57.32	1010.24	76.64	446.48	1	3.082192	0.454822	-0.001897	-0.168222
4	10.82	37.50	1009.23	96.62	473.90	1	-0.481304	-0.345778	-0.001000	0.260699

정규화(Normalize)

# for 문으로 데이터프레임을 뽑아오면 컬럼명이 나오게 됩니다.
# 컬럼명을 키워드로 하여 해당 컬럼의 모든 열을 가져올 수 있습니다.
# 이 부분은 combined[i] 이고, 여기서 i 는 컬럼의 이름 입니다.


# 아래 작업을 모든 컬럼에 대해 반복합니다.
for i in combined_change:
    if i == 'PE':
        continue
    # 한 컬럼에서 최소값과 최대값을 구하고
    # 해당 컬럼의 모든 열 값에서 최소값을 뺍니다.
    # 그리고 이를 최대값으로 나누게 되면, 모든 값들이 0 에서 1 사이의 값으로 정규화 됩니다.
    minimum = combined_change[i].min()
    maximum = combined_change[i].max()
    # 최대값이 아닌, minimum 을 뺀 값으로 나누는 이유는 
    # 모든 열에서 minimum 을 뺀 시점에서 이미 최대값이 minimum 만큼 감소하기 때문입니다.
    combined_change[i] = (combined_change[i] - minimum) / (maximum - minimum)

combined_change.describe()

	AT	V	AP	RH	PE	year	1d_AT	1d_V	1d_AP	1d_RH
count	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000	47840.000000
mean	0.505417	0.515050	0.504060	0.640067	454.365009	0.500000	0.066246	0.265186	0.489677	0.221845
std	0.211109	0.226110	0.146957	0.195706	17.066281	0.353557	0.050898	0.128753	0.125634	0.097108
min	0.000000	0.000000	0.000000	0.000000	420.260000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	0.331445	0.291459	0.401138	0.506267	439.750000	0.250000	0.034131	0.168679	0.405556	0.156678
50%	0.525071	0.475445	0.496164	0.662399	451.550000	0.500000	0.053211	0.244571	0.489235	0.207448
75%	0.677337	0.732740	0.603069	0.794504	468.430000	0.750000	0.081441	0.339624	0.572324	0.268773
max	1.000000	1.000000	1.000000	1.000000	495.760000	1.000000	1.000000	1.000000	1.000000	1.000000

# 각 환경변수들의 한 시간 이전 대비 변화율을 뽑아내었으니
# 이를 모든 변수들과 산점도 행렬을 통해 비교합니다.

# 여기에서 볼 사항은 
# 새롭게 만든 환경변수의 변화율이 어떤 분포를 가지는가 입니다.

from pandas.plotting import *
data = combined_change
scatter(data.drop(columns=['PE']), data['PE'])

# 새롭게 만든 환경변수의 변화율 변수만 따로 뽑아서 보겠습니다.
# 1d_AT 의 경우 우리가 파악하고자 하는 값들이 가장 최소 값 인근으로 쏠려 있음을 확인할 수 있습니다.

features = ['1d_AT', '1d_V', '1d_AP', '1d_RH']
scatter(combined_change[features], combined_change['PE'])

수리 모델링

# 한 쪽으로 쏠려있는 데이터를 고루 퍼지게 만드는 방법에 대해 소개합니다.
# root 를 씌우는 것과 역수를 취하는 방법에 대해 다룹니다.

import numpy as np
target_col = "1d_AT"

combined_change[target_col + "_sqrt"] = np.sqrt(combined_change[target_col] 
                                                - combined_change[target_col].min() + 1)

combined_change[target_col + "_invert_1"] = np.power(combined_change[target_col] 
                                                     - combined_change[target_col].min() + 1, -1)

combined_change[target_col + "_invert_5"] = np.power(combined_change[target_col] 
                                                     - combined_change[target_col].min() + 1, -5)

combined_change[target_col + "_invert_10"] = np.power(combined_change[target_col] 
                                                      - combined_change[target_col].min() + 1, -10)

# -1 승부터 -10 승까지의 데이터를 보면서, 어떤 것을 선택할지 알아봅니다.
# 히스토그램이 정규분포에 가장 가까운 것을 선택합니다.

features = ['1d_AT', '1d_AT_sqrt', '1d_AT_invert_1', '1d_AT_invert_5', '1d_AT_invert_10', 'RH']
scatter(combined_change[features], combined_change['PE'])

# 나머지 파생변수들은 -1승을 취해줍니다.

target_col = "1d_V"
combined_change[target_col + "_invert_1"] = np.power(combined_change[target_col] 
                                                     - combined_change[target_col].min() + 1, -1)

target_col = "1d_RH"
combined_change[target_col + "_invert_1"] = np.power(combined_change[target_col] 
                                                     - combined_change[target_col].min() + 1, -1)

target_col = "1d_AP"
combined_change[target_col + "_invert_1"] = np.power(combined_change[target_col] 
                                                     - combined_change[target_col].min() + 1, -1)

# 원래 있던 변수와 새로 만든 파생변수들을 넣고, 산점도 행렬을 알아봅니다.

features = ['V',  
            'AP',
            'AT', 
            'RH',
            '1d_V_invert_1', 
            '1d_AP_invert_1', 
            '1d_AT_invert_10', 
            '1d_RH_invert_1']

scatter(combined_change[features], combined_change['PE'])

# 원래 있던 변수와, 새롭게 만든 모든 파생변수들을 넣고
# 발전소의 발전량을 알아보는 예측 모델을 만들기

features = ['V',  
            'AP',
            'AT', 
            'RH',
            '1d_V',
            '1d_AP',
            '1d_AT',
            '1d_RH',
            '1d_V_invert_1', 
            '1d_AP_invert_1', 
            '1d_AT_invert_10', 
            '1d_RH_invert_1']

# 데이터를 1000 개를 기점으로 잘라서 학습 / 예측 합니다.
test_slice = 10000

# 슬라이싱을 적용하여 데이터를 특정 기점으로 잘라줍니다.
# 잘랐을 때, 특징 모음의 경우 컬럼의 개수가 동일해야 합니다.

# X_train -> 학습할 데이터의 특징들의 모음
# Y_train -> 학습할 데이터의 레이블

X_train = data[features][:test_slice]
Y_train = data['PE'][:test_slice]
X_train.shape

(10000, 12)

# X_test -> 예측할 데이터의 특징 모음
# Y_test -> 예측할 데이터의 레이블

X_test = data[features][test_slice:]
Y_test = data['PE'][test_slice:]
X_test.shape

(37840, 12)

# 테스트 한 알고리즘을 평가할 점수를 정의합니다.
# RMSE(Root Mean Squared Error) 를 사용합니다.

from sklearn.metrics import make_scorer
import numpy as np

def accuracy(predict, actual):
    predict = np.array(predict)
    actual = np.array(actual)

    difference = actual-predict
    squared = difference ** 2
    root = np.sqrt(squared)
    mean = np.mean(root)

    score = np.mean(difference)

    return score

simple_score = make_scorer(accuracy)
simple_score

make_scorer(accuracy)

# Cross Validation 은 학습한 데이터 안에서의 점수를 평가하는 것입니다.
# model 은 어떤 알고리즘으로 학습할 것인지를 나타냅니다.
# RandomForestRegressor 는 의사결정나무를 여러 개 만들어서 서로 투표하는 방법을 사용합니다.

from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor(random_state=2000, n_estimators=100, n_jobs=8)
cross_val_score(model, 
                X_train, 
                Y_train, 
                cv=20, 
                scoring=simple_score).mean()

 

2.4646049100000007

# model.fit -> 선택한 알고리즘에 학습할 데이터를 넣어서 학습시킨다는 의미입니다.
# model.predict -> 학습한 알고리즘에 예측할 데이터를 넣어서 결과를 예측하는 것입니다.
# accuracy -> 예측한 결과인 predict 와 실제 결과를 비교하여 점수를 환산합니다. 

model.fit(X_train, Y_train)
predict = model.predict(X_test)
accuracy(predict, Y_test)

1.8244322093023257

# 학습(fit) 한 이후에는, 이 예측모델이 어떤 기준에 따라 판별했는지를 알 수 있습니다.
# 이는 의사결정나무 기반 알고리즘의 특징이며, 변수들의 우선순위를 파악하는 데 도움이 됩니다.

print ("Features sorted by their score:")
for i in sorted(zip(map(lambda x: round(x, 4), model.feature_importances_), features), reverse=True):
    print(i)

Features sorted by their score:
(0.8965, 'AT')
(0.0583, 'V')
(0.0143, 'AP')
(0.0124, 'RH')
(0.0028, '1d_RH')
(0.0027, '1d_RH_invert_1')
(0.0025, '1d_AP_invert_1')
(0.0024, '1d_AP')
(0.002, '1d_V_invert_1')
(0.002, '1d_V')
(0.002, '1d_AT_invert_10')
(0.002, '1d_AT')

# 1d_RH 가 1d_RH_invert_1 보다 높은 점수를 얻었으므로, 둘 중에 이를 선택합니다.
# 나머지 파생변수들은 invert 한 것이 그렇지 않은 경우보다 높은 점수를 얻었으므로, 이를 선택합니다.

feature_selected = ['V', 
                    'AP', 
                    'AT', 
                    'RH',
                    '1d_V_invert_1', 
                    '1d_AP_invert_1', 
                    '1d_AT_invert_10', 
                    '1d_RH']

# 어떤 컬럼을 선택할지 정했으므로, 선별된 특징들로 학습 및 검증 데이터를 만들어줍니다.
X_train = data[feature_selected][:test_slice]
Y_train = data['PE'][:test_slice]
X_train.shape

(10000, 8)

X_test = data[feature_selected][test_slice:]
Y_test = data['PE'][test_slice:]
X_test.shape

(37840, 8)

from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor(random_state=2000, n_estimators=100, n_jobs=8)
cross_val_score(model, 
                X_train, 
                Y_train, 
                cv=20, 
                scoring=simple_score).mean()

2.45433224

# 학습할 변수들을 선별하니 점수가 더 좋아졌음을 알 수 있습니다.

model.fit(X_train, Y_train)

predict = model.predict(X_test)
accuracy(predict, Y_test)

1.799808337737842

print ("Features sorted by their score:")
for i in sorted(zip(map(lambda x: round(x, 4), model.feature_importances_), feature_selected), reverse=True):
    print(i)

Features sorted by their score:
(0.8968, 'AT')
(0.0585, 'V')
(0.0146, 'AP')
(0.0128, 'RH')
(0.0052, '1d_RH')
(0.0046, '1d_AP_invert_1')
(0.0038, '1d_AT_invert_10')
(0.0037, '1d_V_invert_1')