[코드 분석] 2023 전력사용량 예측 AI 경진대회

[데이콘] 전력수요예측, XGBoost

j_sean팀 | Private 1위(5.0293) | XGBoost 단일 모형

 

데이터 분석 및 머신러닝을 위한 여러 라이브러리를 import하는 코드

import sys 
import sktime # 시계열 데이터 분석을 위한 Python 패키지
import tqdm as tq # 반복문의 진행 상황을 시각화해주는 라이브러리
import xgboost as xgb # 그래디언트 부스팅 알고리즘을 구현한 라이브러리
# 트리 기반의 앙상블 학습 알고리즘 중 하나
import matplotlib # 데이터 시각화를 위한 파이썬 라이브러리
import seaborn as sns # matplotlib을 기반으로 하는 데이터 시각화 라이브러리
import sklearn as skl # 간단하고 효율적인 데이터 분석과 머신러닝을 위한 라이브러리
import pandas as pd # 데이터 조작 및 분석을 위한 라이브러리
import numpy as np # 수치 연산을 위한 파이썬 라이브러리

 

0. load the libararies

import matplotlib.pyplot as plt # 데이터를 시각화하기 위한 matplotlib의 서브 모듈
from tqdm import tqdm
from sktime.forecasting.model_selection import temporal_train_test_split 
# 시계열 데이터를 train set과 test set으로 분할 ( 시계열 데이터의 시간 정보를 고려 )
from sktime.utils.plotting import plot_series # 시계열 데이터를 그래프로 그리
from xgboost import XGBRegressor # 그래디언트 부스팅 기반의 머신러닝 알고리즘 (회귀 전용)

 

1. preprocessing the data

time series를 일반 regression 문제로 변환하기 위해 시간 관련 변수 추가(월 / 주 / 요일)

전력소비량의 (건물별, 요일별, 시간대별) 평균 / (건물별, 시간대별) 평균 / (건물별, 시간대별, 표준편차) 변수 추가

• 건물별, 요일별, 시간대별 표준편차 / 건물별, 평균 등 여러 통계량 생성 후 몇 개 건물에 테스트, 최종적으로 성능 향상에 도움이 된 위 3개 변수만 추가

공휴일 변수 추가

시간(hour)은 cyclical encoding하여 변수 추가(sin time & cos time) 후 삭제

CDH(Cooling Degree Hour) & THI(불쾌지수) 변수 추가

건물별 모델 생성 시 무의미한 태양광 발전 시설 / 냉방시설 변수 삭제

# 데이터 불러오기
train = pd.read_csv('./data/train.csv', encoding = 'cp949')

## 변수들을 영문명으로 변경
cols = ['num', 'date_time', 'power', 'temp', 'wind','hum' ,'prec', 'sun', 'non_elec', 'solar']
train.columns = cols

# 시간 관련 변수들 생성
date = pd.to_datetime(train.date_time)
train['hour'] = date.dt.hour
train['day'] = date.dt.weekday
train['month'] = date.dt.month
train['week'] = date.dt.weekofyear

## 건물별, 요일별, 시간별 발전량 평균 넣어주기
power_mean = pd.pivot_table(train, values = 'power', index = ['num', 'hour', 'day'], aggfunc = np.mean).reset_index()
tqdm.pandas()
train['day_hour_mean'] = train.progress_apply(lambda x : power_mean.loc[(power_mean.num == x['num']) & (power_mean.hour == x['hour']) & (power_mean.day == x['day']) ,'power'].values[0], axis = 1)

## 건물별 시간별 발전량 평균 넣어주기
power_hour_mean = pd.pivot_table(train, values = 'power', index = ['num', 'hour'], aggfunc = np.mean).reset_index()
tqdm.pandas()
train['hour_mean'] = train.progress_apply(lambda x : power_hour_mean.loc[(power_hour_mean.num == x['num']) & (power_hour_mean.hour == x['hour']) ,'power'].values[0], axis = 1)

## 건물별 시간별 발전량 표준편차 넣어주기
power_hour_std = pd.pivot_table(train, values = 'power', index = ['num', 'hour'], aggfunc = np.std).reset_index()
tqdm.pandas()
train['hour_std'] = train.progress_apply(lambda x : power_hour_std.loc[(power_hour_std.num == x['num']) & (power_hour_std.hour == x['hour']) ,'power'].values[0], axis = 1)

### 공휴일 변수 추가
train['holiday'] = train.apply(lambda x : 0 if x['day']<5 else 1, axis = 1)
train.loc[('2020-08-17'<=train.date_time)&(train.date_time<'2020-08-18'), 'holiday'] = 1

## <https://dacon.io/competitions/official/235680/codeshare/2366?page=1&dtype=recent>
# 시간의 연속성을 나타내기 위한 코드
train['sin_time'] = np.sin(2*np.pi*train.hour/24)
train['cos_time'] = np.cos(2*np.pi*train.hour/24)

## <https://dacon.io/competitions/official/235736/codeshare/2743?page=1&dtype=recent>
train['THI'] = 9/5*train['temp'] - 0.55*(1-train['hum']/100)*(9/5*train['hum']-26)+32

# 각 건물에 대한 온도 데이터를 기반으로 Cold Degree Hour를 계산하고
# 이를 새로운 변수로 추가하여 난방 필요도를 나타내는데 사용한다.
def CDH(xs):
    ys = []
    for i in range(len(xs)):
        if i < 11:
            ys.append(np.sum(xs[:(i+1)]-26))
        else:
            ys.append(np.sum(xs[(i-11):(i+1)]-26))
    return np.array(ys)

cdhs = np.array([])
for num in range(1,61,1):
    temp = train[train['num'] == num]
    cdh = CDH(temp['temp'].values)
    cdhs = np.concatenate([cdhs, cdh])
train['CDH'] = cdhs

train.drop(['non_elec','solar','hour'], axis = 1, inplace = True)
train.head()

 

2. Model : XGBoost

모델은 시계열 데이터에 좋은 성능을 보이는 XGBoost를 선정

# Define SMAPE loss function
def SMAPE(true, pred):
    return np.mean((np.abs(true-pred))/(np.abs(true) + np.abs(pred))) * 100

새로운 objective function은 MSE 대신 SMAPE를 사용하여 모델을 훈련한다.

일반적인 경우, MSE를 objective function으로 사용하면 일부 건물의 예측이 과소추정될 수 있는데, 이로 인해 SMAPE 점수가 높아질 수 있다. 따라서, residual(잔차)이 0보다 클 때, 즉 실제 값보다 낮게 추정될 때 해당 건물의 가중치를 높여 과소 추정되는 건물에 대한 페널티를 높일 수 있다.

 

모델 학습

xgb_params = pd.read_csv('./parameters/hyperparameter_xgb.csv')

xgb_reg = XGBRegressor(n_estimators = 10000, eta = xgb_params.iloc[47,1], min_child_weight = xgb_params.iloc[47,2], 
                       max_depth = xgb_params.iloc[47,3], colsample_bytree = xgb_params.iloc[47,4], 
                       subsample = xgb_params.iloc[47,5], seed=0)

xgb_reg.fit(x_train, y_train, eval_set=[(x_train, y_train), (x_valid, y_valid)],
        early_stopping_rounds=300,
       verbose=False)

 

모델 예측

pred = xgb_reg.predict(x_valid)
pred = pd.Series(pred)
pred.index = np.arange(y_valid.index[0], y_valid.index[-1]+1)
plot_series(y_train, y_valid, pd.Series(pred), markers=[',' , ',', ','])

print('best iterations: {}'.format(xgb_reg.best_iteration))
print('SMAPE : {}'.format(SMAPE(y_valid, pred)))

3. model tuning

다음과 같은 방법으로 하이퍼 파라미터를 튜닝했다.

1. sklearn의 gridsearchCV를 활용해 튜닝
2. XGBoost의 early stopping 기능으로 n_estimators를 튜닝
3. weighted_mse의 alpha값을 튜닝

alpha_list = []
smape_list = []
for i in tqdm(range(60)):
    y = train.loc[train.num == i+1, 'power']
    x = train.loc[train.num == i+1, ].iloc[:, 3:]
    y_train, y_test, x_train, x_test = temporal_train_test_split(y = y, X = x, test_size = 168)
    xgb = XGBRegressor(seed = 0,
                      n_estimators = best_it[i], eta = 0.01, min_child_weight = xgb_params.iloc[i, 2],
                      max_depth = xgb_params.iloc[i, 3], colsample_bytree = xgb_params.iloc[i, 4], subsample = xgb_params.iloc[i, 5])
    
    xgb.fit(x_train, y_train)
    pred0 = xgb.predict(x_test)
    best_alpha = 0
    score0 = SMAPE(y_test,pred0)
    
    for j in [1, 3, 5, 7, 10, 25, 50, 75, 100]:
        xgb = XGBRegressor(seed = 0,
                      n_estimators = best_it[i], eta = 0.01, min_child_weight = xgb_params.iloc[i, 2],
                      max_depth = xgb_params.iloc[i, 3], colsample_bytree = xgb_params.iloc[i, 4], subsample = xgb_params.iloc[i, 5])
        xgb.set_params(**{'objective' : weighted_mse(j)})
    
        xgb.fit(x_train, y_train)
        pred1 = xgb.predict(x_test)
        score1 = SMAPE(y_test, pred1)
        if score1 < score0:
            best_alpha = j
            score0 = score1
    
    alpha_list.append(best_alpha)
    smape_list.append(score0)
    print("building {} || best score : {} || alpha : {}".format(i+1, score0, best_alpha))

seed ensemble

seed의 영향을 제거하기 위해 6개의 seed(0부터 5)별로 훈련, 예측하여 6개 예측값의 평균을 구했다.

preds = np.array([]) 
for i in tqdm(range(60)):
    
    pred_df = pd.DataFrame()   # 시드별 예측값을 담을 data frame
    
    for seed in [0,1,2,3,4,5]: # 각 시드별 예측
        y_train = train.loc[train.num == i+1, 'power']
        x_train, x_test = train.loc[train.num == i+1, ].iloc[:, 3:], test.loc[test.num == i+1, ].iloc[:,1:]
        x_test = x_test[x_train.columns]
        
        xgb = XGBRegressor(seed = seed, n_estimators = best_it[i], eta = 0.01, 
                           min_child_weight = xgb_params.iloc[i, 2], max_depth = xgb_params.iloc[i, 3], 
                           colsample_bytree=xgb_params.iloc[i, 4], subsample=xgb_params.iloc[i, 5])
    
        if xgb_params.iloc[i,6] != 0:  # 만약 alpha가 0이 아니면 weighted_mse 사용
            xgb.set_params(**{'objective':weighted_mse(xgb_params.iloc[i,6])})
        
        xgb.fit(x_train, y_train)
        y_pred = xgb.predict(x_test)
        pred_df.loc[:,seed] = y_pred   # 각 시드별 예측 담기
        
    pred = pred_df.mean(axis=1) # (i+1)번째 건물의 예측 =  (i+1)번째 건물의 각 시드별 예측 평균값
    preds = np.append(preds, pred) 
    
preds = pd.Series(preds)

fig, ax = plt.subplots(60, 1, figsize=(100,200), sharex = True)
ax = ax.flatten()
for i in range(60):
    train_y = train.loc[train.num == i+1, 'power'].reset_index(drop = True)
    test_y = preds[i*168:(i+1)*168]
    ax[i].scatter(np.arange(2040) , train.loc[train.num == i+1, 'power'])
    ax[i].scatter(np.arange(2040, 2040+168) , test_y)
    ax[i].tick_params(axis='both', which='major', labelsize=6)
    ax[i].tick_params(axis='both', which='minor', labelsize=4)
#plt.savefig('./predict_xgb.png')
plt.show()

 

왜 Seed 값을 평균을 냈는지 의문이 들어 GPT에게 물어봤다. 

seed 값은 모델 훈련 과정 중 발생하는 여러 무작위 프로세스의 초기화에 사용된다. 이러한 무작위 프로세스에는 데이터 셔플링, 특성 선택, 분할 포인트의 결정 등이 포함된다.

1. 재현성(Reproducibility): 동일한 seed 값으로 모델을 여러 번 훈련하면, 매번 동일한 결과를 얻을 수 있다. 이는 실험의 재현성을 보장하고, 다른 설정이나 조건을 변경했을 때 결과의 변화를 정확하게 평가하는 데 중요하다.
2. 모델 다양성(Model Diversity): 앙상블 학습에서 여러 모델을 훈련할 때 각각 다른 seed 값을 사용하면, 모델 간의 다양성이 증가합니다. 이 다양성은 모델이 데이터의 다른 측면을 학습하게 하여, 전체 앙상블의 성능을 향상시킬 수 있다. 각 모델이 서로 다른 패턴을 감지하고, 이를 종합함으로써 예측의 정확도와 견고성이 높아진다.

seed 값을 다르게 설정함으로써 생성된 모델들은 초기 가중치 할당, 데이터 샘플링 방식, 트리의 분할 지점 선택 등에서 미묘하게 다른 경로를 따른다. 이러한 차이는 각 모델이 조금씩 다른 예측을 생성하게 만든다. 이들 예측을 평균내거나 다른 방법으로 종합함으로써, 개별 모델의 예측에 포함될 수 있는 무작위 오차나 편향을 상쇄시키는 효과가 있다.

 

나름 합리적인 방법인 것 같아 나중에 활용해보려 한다.

4. post processing

weighted mse와 같은 맥락에서, 과도한 underestimate를 막기 위해 예측값 후처리 진행 예측주로 부터 직전 4주(train set 마지막 28일)의 건물별, 요일별, 시간대별 전력소비량의 최솟값을 구한 뒤, test set의 같은 건물 요일 시간대의 예측값과 비교하여 만약 1번의 최솟값보다 예측값이 작다면 최솟값으로 예측값을 대체해준다. (public score 0.01 , private score 0.08 정도의 성능 향상)

 

---

오차 지표와 랜덤이라는 특성에 대한 깊은 고민을 한 코드라는 점이 느껴졌다. 다른 코드를 읽었을땐 데이터 전처리에 대한 고민한 흔적이 많았었는데 아무래도 초반에 제출된 프로젝트다 보니 결과 값에 좀더 비중을 둔 느낌이다. 다음엔 코드 전처리에 대한 깊은 고찰에 관한 코드를 읽어 봐야겠다.