时间序列特征构造：以电力负荷预测为例讲解（python语言）

个人电气博文目录传送门
学好电气全靠它，个人电气博文目录（持续更新中…）

时间序列特征构造

时间序列问题，首先不管是回归问题，还是分类问题。

一个模型的好坏，决定因素由数据集的大小，特征值的选取和处理，算法。
其中最重要的是特征值的选取和处理。
今天余总来讲解下时间序列的特征构造问题。
该特征构造部分可以用于其他数值数据。

时间序列特征构造分类为：时间特征，时间历史特征，时间交叉特征

时间特征

连续时间：
持续时间，间隔时间

离散时间：
年，季度，季节，月，星期，日，等
节假日，节假日第几天
上午，早上，中午，晚上
年初，年末，月初，月末，周内，周末
是否高峰时段，是否上班

时间历史特征

统计值：
四分位数，中位数，平均值，偏度，峰度，离散系数

同期值：
如今天的8点和昨天的8点，
这周一和上周一

时间交叉特征

类别特征与类别特征：笛卡尔积

连续特征与类别特征：离散笛卡尔积，聚合特征（聚合指聚类后的类别特征）

连续特征与连续特征：一阶差分，二阶差分（即把数据和前后做差值，差值作为特征，一次差值为一阶差分）

时间特征构造代码
使用的包 pandas
个人pandas教程总结
pandas常用函数，个人常用的

时间模块讲解链接
python 日期和时间处理（time，datetime模块讲解）

时间特征构造
原始数据total_df.csv

给出的数据只有时间和负荷，
我们要从时间里提取出年，月，日，节假日，工作日等时间特征。

首先把文档里的时间数据格式由字符串转换到时间

total_df.loc[:, 'record_date'] = pd.to_datetime(total_df['record_date'],format='%Y-%m-%dT%H:%M:%S')

利用python里的time 函数模块，构造强时间指代特征：星期几（dow）、几号（dom）、几月（month）、哪一年（year）。代码如下：


# 几时# 星期几，其中dow为构造的新列
total_df.loc[:, 'dow'] = total_df['record_date'].apply(lambda x: x.dayofweek)# 几号，dom为构造的新列
total_df.loc[:, 'dom'] = total_df['record_date'].apply(lambda x: x.day)# 几月，month为构造的新列
total_df.loc[:, 'month'] = total_df['record_date'].apply(lambda x: x.month)# 几年,year为构造的新列
total_df.loc[:, 'year'] = total_df['record_date'].apply(lambda x: x.year)

利用time函数模块添加周末特征：通过0,1化添加周末特征，周末的特征为1，将文字数值化传入模型，具体代码如下：


#weekend，weekend_sat，weekend_sun 为构造的新列
total_df.loc[:,'weekend'] = 0
total_df.loc[:,'weekend_sat'] = 0
total_df.loc[:,'weekend_sun'] = 0
total_df.loc[(total_df['dow']>4), 'weekend'] = 1
total_df.loc[(total_df['dow']==5), 'weekend_sat'] = 1
total_df.loc[(total_df['dow']==6), 'weekend_sun'] = 1

添加特征：添加上下中旬，其特征分别表示为1/2/3，代码如下：

def period_of_month(day):if day in range(1, 11):return 1if day in range(11, 21):return 2else:return 3
#period_of_month为构造的新列total_df.loc[:, 'period_of_month'] = total_df['dom'].apply(lambda x: period_of_month(x))

添加上半月、下半月时间特征，其特征分别表示为1/2，代码如下：

def period2_of_month(day):if day in range(1, 16):return 1else:return 2
#period2_of_month为构造的新列total_df.loc[:, 'period2_of_month'] = total_df['dom'].apply(lambda x: period2_of_month(x))

利用python 里的time函数模块再自定义函数，添加国庆节假日特征，如果满足月份为10，日期小于8，则为国庆节。代码如下：

total_df.loc[:, 'festival'] = 0total_df.loc[(total_df.month == 10) & (total_df.dom < 8), 'festival'] = 1

天气特征构造
天气啥的一般是文字，（如晴，多云…）文字程序无法直接输入，需要先One-hot化。
one-hot方法有很多，介绍下pandas 包的吧
介绍pd.get_dummies
示例

import pandas as pd
df = pd.DataFrame([['green' , 'A'],['red'   , 'B'],['blue'  , 'A']])df.columns = ['color',  'class']
print('处理前')
print(df)
df1=pd.get_dummies(df)
print('处理后')
print(df1)

一阶差分，二阶差分啥的,pandas shift,diff就能搞定

data['use1']=data['use'].shift(1)#向下平移1
data['use2']=data['use'].shift(2)#向下平移2
data['use3']=data['use'].shift(3)##向下平移3
data['use1差值']=data['use'].diff(1)#前后 t和t-1差值
data['use2差值']=data['use'].diff(2)##t和t-2差值
data= data.dropna()#去空值

统计特征啥的，pandas 有函数了调，不讲解，如果这都不会，自己拿豆腐撞墙吧。

聚合特征：首先把各个样本聚类，然后把类别信息作为特征之一，传入到输入中，作为输入的一部分。

例子：负荷预测

处理好后的数据样子截图

数据已经对天气，节假日进行了处理。

误差函数


#statistic.py
import numpy as np
import math
import scipy.stats as stats#计算预测值真实值的均方误差
#actual 实际值
def meanSquareError(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0total = 0.0for x in range(len(actual)):total += math.pow(actual[x]-pred[x],2)return total/len(actual)
# actual values实际值
#计算mse误差
def mse(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0total = 0.0for x in range(len(actual)):total += math.pow(actual[x]-pred[x],2)return total/(len(actual)*1000000)# 计算预测值与实际值之间的归一化均方根误差(NRMSE)
#
def normRmse(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0sumSquares = 0.0maxY = actual[0]minY = actual[0]for x in range(len(actual)):sumSquares += math.pow(pred[x]-actual[x],2.0)maxY = max(maxY,actual[x])minY = min(minY,actual[x])return math.sqrt(sumSquares/len(actual))/(maxY-minY)# 根据实际值计算预测值的均方根误差(RMSE)def Rmse(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0sumSquares = 0.0for x in range(len(actual)):sumSquares += math.pow(pred[x]-actual[x],2.0)return math.sqrt(sumSquares/len(actual))#计算相对于实际值的预测值的绝对百分误差(MAPE)
def mape(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0total = 0.0for x in range(len(actual)):total += abs((actual[x]-pred[x])/actual[x])return total/len(actual)#计算相对于实际值的预测值的绝对百分误差(MAPE)
def mae(actual,pred):if (not len(actual) == len(pred) or len(actual) == 0):return -1.0total = 0.0for x in range(len(actual)):total += abs(actual[x]-pred[x])return total/len(actual)

预测代码

import math
from tools import statistics#引入误差函数
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
# 将时间序列转换为监督学习问题
#n_in=24时即为将24小时的所有数据去预测24小时后的数据。
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):n_vars = 1 if type(data) is list else data.shape[1]df = DataFrame(data)cols, names = list(), list()# input sequence (t-n, ... t-1)for i in range(n_in, 0, -1):cols.append(df.shift(i))#当names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]# forecast sequence (t, t+1, ... t+n)for i in range(0, n_out):cols.append(df.shift(-i))#向上平移即将下一天当前时间的actual功率平移当当前天，作为输出标签if i == 0:names += [('var%d(t)' % (j+1)) for j in range(n_vars)]else:names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]# put it all togetheragg = concat(cols, axis=1)agg.columns = names# drop rows with NaN valuesif dropnan:agg.dropna(inplace=True)return agg
# 加载数据
dataset = read_csv('load.csv', header=0, index_col=0)
values = dataset.values
# 整数编码
# encoder = LabelEncoder()
# values[:,4] = encoder.fit_transform(values[:,4])
# ensure all data is float
values = values.astype('float32')
# normalize features 归一化特征
# scaler = MinMaxScaler(feature_range=(0, 1))
# scaled = scaler.fit_transform(values)
matrix_load = np.array(values)
matrix_load[:,0] = matrix_load[:,0]/1000#daily load 加载日负载数据，并将第一列功率每一个数除以1000# 异常值处理
k = 0
for j in range(0, matrix_load.shape[0]):
##如何一个数与前一个相差2（即2000）和与后一个相差2（2000）if(abs(matrix_load[j,0]-matrix_load[j-1,0])>2 and abs(matrix_load[j,0]-matrix_load[j+1,0])>2):k = k + 1##则这个数=（前一个+后一个）/2+(前一天的前一个+前一天的后一个)/2matrix_load[j,0] = (matrix_load[j - 1,0] + matrix_load[j + 1,0]) / 2 + matrix_load[j - 24,0] - matrix_load[j - 24 - 1,0] / 2sum = 0num = 0for t in range(1,8):if(j - 24*t >= 0):num = num + 1sum = sum + matrix_load[j - 24*t,0]#前一周以内每日同时数据之和（如星期一的4点，星期二的4点...)if((j + 24*t) < matrix_load.shape[0]):num = num + 1sum = sum + matrix_load[j + 24*t,0]#前一周的+后一周每日同时数据之和sum = sum / num#前后一周每日同时数据平均值# #如果当前数据与平均值相差超过3，则当前数据=平均值加或减3if(abs(matrix_load[j,0] - sum)>3):k = k + 1if(matrix_load[j,0] > sum): matrix_load[j,0] = sum + 3else: matrix_load[j,0] = sum - 3
print(k)
# print(matrix_load[1,0])
# shift all data by mean 去均值
shifted_value = matrix_load[:,0].mean()
matrix_load[:,0] -= shifted_value
for i in range(1, 9):matrix_load[:, i] = matrix_load[:, i] / 10shifted_valuei = matrix_load[:,i].mean()matrix_load[:,i] -= shifted_valuei
print(matrix_load.shape)
# frame as supervised learning
before = 24
end = 1
reframed = series_to_supervised(matrix_load, before, end)
print(reframed.shape)
# drop columns we don't want to predict
for i in range(76):reframed.drop(reframed.columns[[-1]], axis=1, inplace=True)
print(reframed.shape)
print(reframed.head())
# 划分训练集测试集
values = reframed.values
n_train_hours = values.shape[0] - 166*24
train = values[:n_train_hours, :]
test = values[n_train_hours:, :]
print('train:',n_train_hours,'test:',166*24)
# 划分输入和输出
train_X, train_y = train[:, :-end], train[:, -end]
print(train_X.shape, train_y.shape)
test_X, test_y = test[:, :-end], test[:, -end]
#改变数据形状 3D [samples, timesteps, features]
train_X = train_X.reshape((train_X.shape[0], before*77))
test_X = test_X.reshape((test_X.shape[0], before*77))
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)
# svr
# kernel='linear'时，为线性核，C越大分类效果越好，但有可能会过拟合（defaul C=1）
# kernel='rbf'时（default），为高斯核radial basis function，gamma值越小，分类界面越连续；gamma值越大，分类界面越“散”，分类效果越好，但有可能会过拟合
kernelList = ["rbf"]
names = ["true","radial basis"]
preds = []
preds.append(test_y)
for i in range(len(kernelList)):clf = svm.SVR(C=2.0, kernel=kernelList[i])clf.fit(train_X, train_y)predicted_values = clf.predict(test_X)mape = statistics.mape((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)print('MAPE is ', mape)mae = statistics.mae((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)print('MAE is ', mae)mse = statistics.meanSquareError((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)print('MSE is ', mse)rmse = math.sqrt(mse)print('RMSE is ', rmse)nrmse = statistics.normRmse((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)print('NRMSE is ', nrmse)preds.append(predicted_values)
# 结果
fig = plt.figure()
colors = ["g","r","b","c","m","y","k","w"]
legendVars = []
for j in range(len(preds)):print(j)x, = plt.plot(preds[j]+shifted_value, color=colors[j])legendVars.append(x)
plt.title("Svm (rbf) Multi")
plt.xlabel('Time（hour）')
plt.ylabel('Electricity load (*1e3 kW)')
plt.legend(legendVars, names)
plt.show()
fig.savefig('svr_mul_result.jpg', bbox_inches='tight')