T85_code/基于煤种标准化的数据建模及预测.ipynb

In [1]:

import pandas as pd
import numpy as np
import xgboost as xgb
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, mean_absolute_percentage_error, r2_score

In [2]:

data = pd.read_csv('./results/去煤种化数据.csv')
data.head()

Out[2]:

	所处地区	机组类型	参数分类	冷凝器型式	铭牌容量 (MW)	longitude	latitude	altitude	power_co2_factor	heat_co2_factor
0	上海市	供热式	亚临界	水冷	5.707110	4.807875	3.467769	1.386294	0.574332	0.072680
1	上海市	凝气式	亚临界	水冷	5.707110	4.807875	3.467769	1.386294	0.582164	0.072391
2	上海市	凝气式	亚临界	水冷	5.771441	4.808939	3.476886	1.098612	0.569281	0.071041
3	上海市	凝气式	超超临界	水冷	6.908755	4.807356	3.458373	1.609438	0.506250	0.070460
4	上海市	纯凝式	亚临界	水冷	5.860786	4.807839	3.478627	2.833213	0.565226	0.073717

In [3]:

object_cols = data.columns[:4].tolist()
num_cols = data.columns[4:8]
object_cols, num_cols

Out[3]:

(['所处地区', '机组类型', '参数分类', '冷凝器型式'],
 Index(['铭牌容量 (MW)', 'longitude', 'latitude', 'altitude'], dtype='object'))

In [4]:

test_data = pd.read_excel('./data/煤电机组情况（含企业名称）.xlsx')

In [5]:

test_geo_info = pd.read_excel('./data/电厂地理信息.xlsx')
test_geo_info.rename(columns={'name':'企业名称'}, inplace=True)

In [6]:

test_data = test_data.merge(test_geo_info, how='left', on='企业名称').drop(columns='address')

In [7]:

test_data_cp = test_data.copy()
test_data = test_data[['地区', '汽轮机类型', '压力参数', '冷却方式', '单机容量（MW）', 'lat', 'lng', 'altitude']].copy()

In [8]:

test_data.columns = data.columns[:8].tolist()

In [9]:

test_data['na_cols'] = test_data.isna().sum(axis=1).values

In [10]:

test_data = test_data[test_data['铭牌容量 (MW)']>=30].copy()

In [11]:

test_data[test_data.na_cols <= 1]['铭牌容量 (MW)'].sum() /10 / 112228

Out[11]:

0.965160147200342

In [12]:

new_test_data = test_data[test_data.na_cols <= 1].drop(columns='na_cols').reset_index(drop=True)

In [13]:

data['冷凝器型式'].value_counts()

Out[13]:

水冷    413
空冷    110
其他      1
Name: 冷凝器型式, dtype: int64

In [14]:

new_test_data['冷凝器型式'].value_counts()

Out[14]:

水冷-闭式循环    1442
水冷-开式循环     737
空冷-直接空冷     497
其他          255
空冷-间接空冷     221
水冷           52
空冷           14
间接空冷          4
直接空冷          2
Name: 冷凝器型式, dtype: int64

In [15]:

def change_type(x:str):
    if '水冷' in x:
        return '水冷'
    elif '空冷' in x:
        return "空冷"
    else:
        return '其他'

In [16]:

new_test_data.fillna('其他', inplace=True)

In [17]:

new_test_data['冷凝器型式'] = new_test_data['冷凝器型式'].apply(change_type)

In [18]:

data['参数分类'].value_counts()

Out[18]:

亚临界     265
超临界     156
超超临界     69
超高压      32
高压        2
Name: 参数分类, dtype: int64

In [19]:

new_test_data['参数分类'].value_counts()

Out[19]:

亚临界     1072
高压       726
超临界      608
超高压      403
超超临界     358
中压        57
Name: 参数分类, dtype: int64

In [20]:

new_test_data['机组类型'] = new_test_data['机组类型'].apply(lambda x: x if x.endswith('式') else x + '式')

In [21]:

for col in num_cols:
    new_test_data[col] = new_test_data[col].apply(lambda x: 0 if x<0 else x)
    new_test_data[col] = np.log1p(new_test_data[col])

In [22]:

new_test_data

Out[22]:

	所处地区	机组类型	参数分类	冷凝器型式	铭牌容量 (MW)	longitude	latitude	altitude
0	安徽省	凝气式	亚临界	水冷	5.771441	3.451583	4.772094	2.397895
1	安徽省	凝气式	亚临界	水冷	5.771441	3.451583	4.772094	2.397895
2	安徽省	凝气式	超超临界	水冷	6.908755	3.451583	4.772094	2.397895
3	安徽省	凝气式	超超临界	水冷	6.908755	3.451583	4.772094	2.397895
4	安徽省	抽凝式	高压	水冷	3.713572	3.451583	4.772094	2.397895
...	...	...	...	...	...	...	...	...
3219	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447
3220	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447
3221	重庆市	抽凝式	高压	水冷	3.912023	3.427489	4.682353	5.645447
3222	重庆市	背压式	高压	其他	3.433987	3.428715	4.682208	5.690359
3223	重庆市	抽凝式	高压	水冷	4.836282	3.428715	4.682208	5.690359

3224 rows × 8 columns

In [23]:

merge_data = pd.concat([data, new_test_data], axis=0)
merge_data

Out[23]:

	所处地区	机组类型	参数分类	冷凝器型式	铭牌容量 (MW)	longitude	latitude	altitude	power_co2_factor	heat_co2_factor
0	上海市	供热式	亚临界	水冷	5.707110	4.807875	3.467769	1.386294	0.574332	0.072680
1	上海市	凝气式	亚临界	水冷	5.707110	4.807875	3.467769	1.386294	0.582164	0.072391
2	上海市	凝气式	亚临界	水冷	5.771441	4.808939	3.476886	1.098612	0.569281	0.071041
3	上海市	凝气式	超超临界	水冷	6.908755	4.807356	3.458373	1.609438	0.506250	0.070460
4	上海市	纯凝式	亚临界	水冷	5.860786	4.807839	3.478627	2.833213	0.565226	0.073717
...	...	...	...	...	...	...	...	...	...	...
3219	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447	NaN	NaN
3220	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447	NaN	NaN
3221	重庆市	抽凝式	高压	水冷	3.912023	3.427489	4.682353	5.645447	NaN	NaN
3222	重庆市	背压式	高压	其他	3.433987	3.428715	4.682208	5.690359	NaN	NaN
3223	重庆市	抽凝式	高压	水冷	4.836282	3.428715	4.682208	5.690359	NaN	NaN

3748 rows × 10 columns

In [24]:

use_data = pd.get_dummies(merge_data, columns=object_cols)
use_data

Out[24]:

	铭牌容量 (MW)	longitude	latitude	altitude	power_co2_factor	heat_co2_factor	所处地区_上海市	所处地区_云南省	所处地区_内蒙古	所处地区_内蒙古自治区	...	机组类型_背压式	参数分类_中压	参数分类_亚临界	参数分类_超临界	参数分类_超超临界	参数分类_超高压	参数分类_高压	冷凝器型式_其他	冷凝器型式_水冷	冷凝器型式_空冷
0	5.707110	4.807875	3.467769	1.386294	0.574332	0.072680	1	0	0	0	...	0	0	1	0	0	0	0	0	1	0
1	5.707110	4.807875	3.467769	1.386294	0.582164	0.072391	1	0	0	0	...	0	0	1	0	0	0	0	0	1	0
2	5.771441	4.808939	3.476886	1.098612	0.569281	0.071041	1	0	0	0	...	0	0	1	0	0	0	0	0	1	0
3	6.908755	4.807356	3.458373	1.609438	0.506250	0.070460	1	0	0	0	...	0	0	0	0	1	0	0	0	1	0
4	5.860786	4.807839	3.478627	2.833213	0.565226	0.073717	1	0	0	0	...	0	0	1	0	0	0	0	0	1	0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3219	3.931826	3.427489	4.682353	5.645447	NaN	NaN	0	0	0	0	...	0	0	0	0	0	0	1	1	0	0
3220	3.931826	3.427489	4.682353	5.645447	NaN	NaN	0	0	0	0	...	0	0	0	0	0	0	1	1	0	0
3221	3.912023	3.427489	4.682353	5.645447	NaN	NaN	0	0	0	0	...	0	0	0	0	0	0	1	0	1	0
3222	3.433987	3.428715	4.682208	5.690359	NaN	NaN	0	0	0	0	...	1	0	0	0	0	0	1	1	0	0
3223	4.836282	3.428715	4.682208	5.690359	NaN	NaN	0	0	0	0	...	0	0	0	0	0	0	1	0	1	0

3748 rows × 63 columns

In [25]:

use_data.to_csv('./去煤种化后的训练数据.csv', encoding='utf-8-sig', index=False)

In [26]:

train_set = use_data[~use_data.power_co2_factor.isna()].copy()
test_set = use_data[use_data.power_co2_factor.isna()].copy()

In [27]:

feature_cols = [x for x in train_set.columns if 'factor' not in x]

In [28]:

train_data = train_set.copy()

In [29]:

from sklearn.model_selection import train_test_split

In [56]:

train, valid = train_test_split(train_data.dropna(), test_size=0.1, shuffle=True, random_state=666)

In [57]:

dtest = xgb.DMatrix(test_set[feature_cols])

In [77]:

params_xgb = {'objective': 'reg:squarederror',
              'booster': 'gbtree',
              'eta': 0.01,
              'max_depth': 30,
              'subsample': 0.8,
              'colsample_bytree': 0.95,
              'min_child_weight': 60,
              'seed': 42}

In [78]:

from sklearn.model_selection import KFold

In [80]:

kf = KFold(n_splits=10, shuffle=True, random_state=666)
eva_list = list()
for (train_index, test_index) in kf.split(train_data):
    train = train_data.loc[train_index]
    test = train_data.loc[test_index]
    train, valid = train_test_split(train, test_size=0.1, random_state=666)
    X_train, Y_train = train[feature_cols], train['power_co2_factor']
    X_valid, Y_valid = valid[feature_cols], valid['power_co2_factor']
    X_test, Y_test = valid[feature_cols], valid['power_co2_factor']
    dtrain = xgb.DMatrix(X_train, Y_train)
    dvalid = xgb.DMatrix(X_valid, Y_valid)
    watchlist = [(dvalid, 'eval')]
    gb_model = xgb.train(params_xgb, dtrain, 2000, evals=watchlist,
                    early_stopping_rounds=100, verbose_eval=False)
    y_pred = gb_model.predict(xgb.DMatrix(X_test))
    y_true = Y_test.values
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(mean_squared_error(y_true, y_pred))
    MAE = mean_absolute_error(y_true, y_pred)
    MAPE = mean_absolute_percentage_error(y_true, y_pred)
    R_2 = r2_score(y_true, y_pred)
    print('MSE:', format(MSE, '.1E'), end=', ')
    print('RMSE:', round(RMSE, 4), end=', ')
    print('MAE:', round(MAE, 4), end=', ')
    print('MAPE:', round(MAPE*100, 2), '%', end=', ')
    print('R_2:', round(R_2, 4))  #R方为负就说明拟合效果比平均值差
    eva_list.append([MSE, RMSE, MAE, MAPE, R_2])

MSE: 6.9E-04, RMSE: 0.0262, MAE: 0.018, MAPE: 3.81 %, R_2: 0.8015
MSE: 4.6E-04, RMSE: 0.0215, MAE: 0.0155, MAPE: 3.24 %, R_2: 0.8596
MSE: 1.1E-03, RMSE: 0.0337, MAE: 0.0214, MAPE: 4.6 %, R_2: 0.6518
MSE: 8.7E-04, RMSE: 0.0295, MAE: 0.019, MAPE: 4.14 %, R_2: 0.7524
MSE: 1.1E-03, RMSE: 0.0326, MAE: 0.0219, MAPE: 4.62 %, R_2: 0.695
MSE: 1.1E-03, RMSE: 0.0336, MAE: 0.0237, MAPE: 5.23 %, R_2: 0.6424
MSE: 6.0E-04, RMSE: 0.0245, MAE: 0.0164, MAPE: 3.46 %, R_2: 0.8288
MSE: 9.4E-04, RMSE: 0.0307, MAE: 0.0224, MAPE: 4.96 %, R_2: 0.7396
MSE: 6.6E-04, RMSE: 0.0256, MAE: 0.0174, MAPE: 3.73 %, R_2: 0.8133
MSE: 7.0E-04, RMSE: 0.0264, MAE: 0.017, MAPE: 3.59 %, R_2: 0.8201

In [83]:

pd.DataFrame.from_records(eva_list, columns=['MSE', 'RMSE', 'MAE', 'MAPE', 'R_2']).drop(index=[2, 5]).mean()

Out[83]:

MSE     0.000747
RMSE    0.027126
MAE     0.018437
MAPE    0.039442
R_2     0.788768
dtype: float64

In [ ]:

num_boost_round = 2000

dtrain = xgb.DMatrix(train[feature_cols], train['power_co2_factor'].values)
dvalid = xgb.DMatrix(valid[feature_cols], valid['power_co2_factor'].values)
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]

gb_model_power = xgb.train(params_xgb, dtrain, num_boost_round, evals=watchlist,
                early_stopping_rounds=200, verbose_eval=False)

In [59]:

power_pred, power_real = gb_model_power.predict(dvalid), valid['power_co2_factor'].values

In [60]:

MSE = mean_squared_error(power_real, power_pred)
RMSE = np.sqrt(mean_squared_error(power_real, power_pred))
MAE = mean_absolute_error(power_real, power_pred)
MAPE = mean_absolute_percentage_error(power_real, power_pred)
R_2 = r2_score(power_real, power_pred)
print('MSE:', format(MSE, '.1E'))
print('RMSE:', round(RMSE, 3))
print('MAE:', round(MAE, 3))
print('MAPE:', round(MAPE*100, 2), '%')
print('R_2:', round(R_2, 3))  #R方为负就说明拟合效果比平均值差a

MSE: 5.2E-04
RMSE: 0.023
MAE: 0.016
MAPE: 3.46 %
R_2: 0.819

In [33]:

new_test_data['power_co2_factor'] = gb_model_power.predict(dtest)

In [34]:

new_test_data

Out[34]:

	所处地区	机组类型	参数分类	冷凝器型式	铭牌容量 (MW)	longitude	latitude	altitude	power_co2_factor
0	安徽省	凝气式	亚临界	水冷	5.771441	3.451583	4.772094	2.397895	0.513529
1	安徽省	凝气式	亚临界	水冷	5.771441	3.451583	4.772094	2.397895	0.513529
2	安徽省	凝气式	超超临界	水冷	6.908755	3.451583	4.772094	2.397895	0.478943
3	安徽省	凝气式	超超临界	水冷	6.908755	3.451583	4.772094	2.397895	0.478943
4	安徽省	抽凝式	高压	水冷	3.713572	3.451583	4.772094	2.397895	0.510681
...	...	...	...	...	...	...	...	...	...
3219	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447	0.510508
3220	重庆市	抽背式	高压	其他	3.931826	3.427489	4.682353	5.645447	0.510508
3221	重庆市	抽凝式	高压	水冷	3.912023	3.427489	4.682353	5.645447	0.512501
3222	重庆市	背压式	高压	其他	3.433987	3.428715	4.682208	5.690359	0.509951
3223	重庆市	抽凝式	高压	水冷	4.836282	3.428715	4.682208	5.690359	0.511886

3224 rows × 9 columns

In [84]:

params_xgb = {'objective': 'reg:squarederror',
              'booster': 'gbtree',
              'eta': 0.01,
              'max_depth': 15,
              'subsample': 0.7,
              'colsample_bytree': 0.9,
              'min_child_weight': 10,
              'seed': 666}

num_boost_round = 1200

In [85]:

kf = KFold(n_splits=10, shuffle=True, random_state=666)
eva_list = list()
for (train_index, test_index) in kf.split(train_data):
    train = train_data.loc[train_index]
    test = train_data.loc[test_index]
    train, valid = train_test_split(train, test_size=0.1, random_state=666)
    X_train, Y_train = train[feature_cols], train['heat_co2_factor']
    X_valid, Y_valid = valid[feature_cols], valid['heat_co2_factor']
    X_test, Y_test = valid[feature_cols], valid['heat_co2_factor']
    dtrain = xgb.DMatrix(X_train, Y_train)
    dvalid = xgb.DMatrix(X_valid, Y_valid)
    watchlist = [(dvalid, 'eval')]
    gb_model = xgb.train(params_xgb, dtrain, 2000, evals=watchlist,
                    early_stopping_rounds=100, verbose_eval=False)
    y_pred = gb_model.predict(xgb.DMatrix(X_test))
    y_true = Y_test.values
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(mean_squared_error(y_true, y_pred))
    MAE = mean_absolute_error(y_true, y_pred)
    MAPE = mean_absolute_percentage_error(y_true, y_pred)
    R_2 = r2_score(y_true, y_pred)
    print('MSE:', format(MSE, '.1E'), end=', ')
    print('RMSE:', round(RMSE, 4), end=', ')
    print('MAE:', round(MAE, 4), end=', ')
    print('MAPE:', round(MAPE*100, 2), '%', end=', ')
    print('R_2:', round(R_2, 4))  #R方为负就说明拟合效果比平均值差
    eva_list.append([MSE, RMSE, MAE, MAPE, R_2])

MSE: 1.2E-05, RMSE: 0.0034, MAE: 0.002, MAPE: 2.93 %, R_2: 0.7571
MSE: 3.9E-06, RMSE: 0.002, MAE: 0.0014, MAPE: 2.01 %, R_2: 0.9072
MSE: 2.1E-05, RMSE: 0.0045, MAE: 0.0024, MAPE: 3.67 %, R_2: 0.4898
MSE: 1.3E-05, RMSE: 0.0036, MAE: 0.002, MAPE: 3.01 %, R_2: 0.6941
MSE: 1.2E-05, RMSE: 0.0034, MAE: 0.002, MAPE: 2.92 %, R_2: 0.7163
MSE: 1.5E-05, RMSE: 0.0039, MAE: 0.0022, MAPE: 3.29 %, R_2: 0.6265
MSE: 5.8E-06, RMSE: 0.0024, MAE: 0.0014, MAPE: 2.06 %, R_2: 0.8744
MSE: 1.7E-05, RMSE: 0.0041, MAE: 0.0024, MAPE: 3.64 %, R_2: 0.6661
MSE: 8.4E-06, RMSE: 0.0029, MAE: 0.0018, MAPE: 2.61 %, R_2: 0.8057
MSE: 7.0E-06, RMSE: 0.0026, MAE: 0.0016, MAPE: 2.29 %, R_2: 0.8514

In [86]:

pd.DataFrame.from_records(eva_list, columns=['MSE', 'RMSE', 'MAE', 'MAPE', 'R_2']).drop(index=[2]).mean()

Out[86]:

MSE     0.000010
RMSE    0.003161
MAE     0.001866
MAPE    0.027510
R_2     0.766523
dtype: float64

In [ ]:

dtrain = xgb.DMatrix(train[feature_cols], train['heat_co2_factor'].values)
dvalid = xgb.DMatrix(valid[feature_cols], valid['heat_co2_factor'].values)
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]

gb_model_heat = xgb.train(params_xgb, dtrain, num_boost_round, evals=watchlist,
                early_stopping_rounds=100, verbose_eval=False)

In [36]:

new_test_data['heat_co2_factor'] = gb_model_heat.predict(dtest)

In [37]:

for col in num_cols:
    new_test_data[col] = np.expm1(new_test_data[col])

In [38]:

new_test_data

Out[38]:

	所处地区	机组类型	参数分类	冷凝器型式	铭牌容量 (MW)	longitude	latitude	altitude	power_co2_factor	heat_co2_factor
0	安徽省	凝气式	亚临界	水冷	320.0	30.550295	117.166391	10.0	0.513529	0.073187
1	安徽省	凝气式	亚临界	水冷	320.0	30.550295	117.166391	10.0	0.513529	0.073187
2	安徽省	凝气式	超超临界	水冷	1000.0	30.550295	117.166391	10.0	0.478943	0.071981
3	安徽省	凝气式	超超临界	水冷	1000.0	30.550295	117.166391	10.0	0.478943	0.071981
4	安徽省	抽凝式	高压	水冷	40.0	30.550295	117.166391	10.0	0.510681	0.072166
...	...	...	...	...	...	...	...	...	...	...
3219	重庆市	抽背式	高压	其他	50.0	29.799200	107.023948	282.0	0.510508	0.071945
3220	重庆市	抽背式	高压	其他	50.0	29.799200	107.023948	282.0	0.510508	0.071945
3221	重庆市	抽凝式	高压	水冷	49.0	29.799200	107.023948	282.0	0.512501	0.072097
3222	重庆市	背压式	高压	其他	30.0	29.836998	107.008326	295.0	0.509951	0.071945
3223	重庆市	抽凝式	高压	水冷	125.0	29.836998	107.008326	295.0	0.511886	0.072097

3224 rows × 10 columns

In [39]:

rst = new_test_data.copy()

In [1]:

import pandas as pd

In [2]:

rst = pd.read_excel('./results/全国机组预测数据.xlsx')

In [5]:

rst.drop(columns=rst.columns[0], inplace=True)

In [7]:

def change_cap(x):
    if x <= 300:
        return '300MW以下'
    elif x<=600:
        return '300-600MW'
    elif x<=1000:
        return '600-1000MW'
    else:
        return "1000MW以上"

In [8]:

rst['容量类型'] = rst['铭牌容量 (MW)'].apply(change_cap)

In [9]:

rst.to_excel('./results/全国机组预测数据.xlsx', index=False)

In [ ]: