用python绘制负荷曲线 python 负荷预测

这是我之前工作做的一个项目

import os
import pandas as pd
import numpy
path = "E:/工作/负荷预测/历史负荷数据-每天" #文件夹目录
files= os.listdir(path) #得到文件夹下的所有文件名称
data_history = pd.DataFrame()
for file in files: #遍历文件夹
if not os.path.isdir(file): #判断是否是文件夹，不是文件夹才打开
f = open(path+"/"+file)
try:
df = pd.read_csv(f)
except UnicodeDecodeError:
print(path+"/"+file)
data_history = pd.concat([data_history,df],ignore_index=True)
data_history = data_history.sort_values(by=['bs','time'])
import sys
import warnings
from pandas.plotting import register_matplotlib_converters
import numpy as np
from datetime import datetime,timedelta
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.arima_model import ARMA
import statsmodels.api as sm
from sklearn.metrics import mean_squared_error
import matplotlib.pylab as plt
import seaborn as sea
%matplotlib inline
from matplotlib.pylab import rcParams
#绘制图像的背景长高
rcParams['figure.figsize'] = 50,22
#用来正常显示中文标签
rcParams['font.sans-serif']=['SimHei']
#用来正常显示负号
rcParams['axes.unicode_minus'] = False
#所有警告只出现一次，避免警告占用太多输出页面，在进行循环时这项设置尤其有用
warnings.filterwarnings(action='once')
data_history = data_history.pivot(index='bs',columns='time',values='value')
#很多台区在2017-08-14之前都存在负荷空值，要把这些时间段去掉
data_history = data_history.iloc[:,41:]
#将索引变为列，方便之后格式转换
data_time_series.index.names=['time_series']
data_time_series=data_time_series.reset_index()
data_time_series.head()
#将含有时序数据的字段转化为datetime64格式
# dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m-%d')
# dateparse('2017/08/01')
#data_time_series['time_series'] = data_time_series['time_series'].apply(dateparse)
date_index = pd.date_range('14/8/2017', '29/8/2019',freq='D')
#检验我们截取的历史负荷的这个时间段是否有遗漏的日期
print(date_index[~date_index.isin(data_time_series['time_series'])])
#将日期索引用datetime64格式的日期代替（python的日期格式我也不太懂，但是用这种格式的日期没报错，我就用这个格式了）
data_time_series['time_series'] = pd.date_range('14/8/2017', '29/8/2019')
data_time_series.set_index('time_series',inplace=True)
data_time_series.head()原始时间序列图

平稳性检测与平稳化处理

ARMA模型要求输入的时间序列类型数据具有平稳性，所以在进行时间序列预测之前，要对时间序列数据进行平稳性的检验和处理。

ADF是一种常用的单位根检验方法，它的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。

通常的平稳化处理方法主要有：对数变换、移动平均、指数平滑、以及差分。其中差分的平稳化处理效果最好，我们在这里采用差分的方法来平稳化时间序列。（差分过程写主程序里了）

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#进行单位根检验:
print('单位根检验:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
if dfoutput.loc['p-value']>0.01:
print('p值大于0.01，不能拒绝原假设，说明时间序列是非平稳的')
else:
print('p值小于0.01，拒绝原假设，说明时间序列是平稳的')
return dfoutput.loc['p-value']<0.01
#返回时间序列是否平稳的判断值，我们需要平稳的时间序列，平稳的时间序列才有可预测性

'''单位根的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。'''ADF检验，检验结果为平稳

建模并确定模型的最佳阶数

将平稳化的趋势因素和随机因素代入ARMA模型之中，但要做出最终的预测，还需要确定模型的阶数(p, q)。一个好的模型通常要求残差序列方差较小，同时模型页相对简单，即要求阶数较低。因此我们需要一个准则来比较不同阶数的模型之间的优劣，从而确定最合适的阶数。

贝叶斯信息准则：定义

使得BIC达到最小值的p即为该准则下的最优AR模型的阶数。

我们采用一种循环的方式，在1-20的范围内寻找能令模型的bic属性最小的阶数值，则此阶数值就是最合适的模型阶数。

#利用贝叶斯信息准则（BIC）寻找最佳阶数，使得BIC达到最小值的(p, q)即为该准则下的最优模型的阶数
def proper_model(data_ts, maxLag):
init_bic = sys.maxsize
init_p = 0
init_q = 0
init_properModel = None
for p in np.arange(maxLag):
for q in np.arange(maxLag):
model = ARMA(data_ts, order=(p, q),freq=data_ts.index.inferred_freq)
try:#
results_ARMA = model.fit(disp=-1, method='css')
except:
continue#忽略所有异常
bic = results_ARMA.bic
if bic < init_bic:
init_p = p
init_q = q
init_properModel = results_ARMA
init_bic = bic
return [init_bic, init_p, init_q, init_properModel]

预测

由于季节性因素是以年为单位循环变动的，那么对季节性因素的预测只需要简单的将上一年的负荷数据照搬到下一年即可。

假如分解出来的趋势因素和随机因素又做了一次平稳化处理的话，则要将趋势因素、随机因素的预测结果做相应的还原。然后将季节性因素、趋势因素、随机因素的预测结果加到一起，就是完整的预测结果

#样本外预测
def compare_ARIMA_modes(series,order):
history_f = [x for x in series]
series_f = series
model = ARIMA(history_f, order=order)
try:
model_fit = model.fit(disp=-1)
except:
flag = pd.Series()
return flag
yhat_f = model_fit.forecast(steps=156)[0]
#附加新元素时也加上一个单位的索引
for t in range(156,1,-1):
series_f = series_f.append(pd.Series(yhat_f[-t],index=[series_f.index[-1]+timedelta(days=7)]))
return series_f[series.index[-1]:series_f.index[-1]]
#根据增长率反推“S”型曲线的阶段
def SCurveStage(y):
#如果增长率小于等于零，说明这个台区的负荷已经达到峰值，不会再增长了，完成了“S”型曲线的增长
if y <=0:
t = float('inf')
return t
t = - np.log((-y)/(y+1-np.e))
return t
#根据所得到的“S”型曲线的阶段，进行增长率预测
def SCurveForecast(t,series):
history_f = [x for x in series]
series_f = pd.Series()
start_time = series.index[-1]
start = series[-1]
#附加新元素时也加上一个单位的索引
for i in range(1,156):
y = (1 + np.exp(1-t))/(1 + np.exp(1-t-1/52))-1
t = t + (1/52)
start = start*(1+y)
series_f = series_f.append(pd.Series(start,index=[start_time+timedelta(weeks=i)]))
return series_f
#季节变动每一年都是一样的，预测季节因素只需要将前一年的数据复制到新一年即可
import copy
def season_forcast_year(seasonal):
seasonal_f = copy.deepcopy(seasonal[-52:-1])
seasonal_f.index = seasonal_f.index+timedelta(weeks=156)
return seasonal_f
industry_type = pd.read_excel('C:/Users/GHZXK/Desktop/预测/台区行业分类.xlsm',sheet_name='基础行业',index_col=0)
industry_type_bs = list(industry_type.index)
df_forecast = pd.DataFrame()

分解时间序列

电力负荷时间序列数据具有非常鲜明的季节性特征和趋势特征，每一个年度的用电数据在节假日的波动上往往具有很大的相似性，而后一年的电力数据往往要比前一年的电力数据高一些。

季节性因素以年为周期，而趋势因素则有很强的线性特征，这两者相对比较好预测。时间序列数据分离出季节性因素和趋势因素后还剩下随机因素，分解后可以避免三种变动因素之间互相影响，并且可以单独测定每种变动的影响程度，从而提高预测精度。因此，我们采用季节分解法来对电力负荷时间序列数据进行分解。

def forecast_decomposition_plt(tqmc,ts,decomposition):
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.figure(figsize=(50,22))
plt.subplot(411)
plt.plot(ts, label='原始')
# 生成网格
plt.grid()
plt.legend(loc='best',fontsize = 60)
plt.savefig(tqmc+'--原始历史负荷.png',bbox_inches='tight')
plt.subplot(412)
#让趋势的纵坐标与原始时间序列一致，避免放大趋势的波动
plt.ylim((min(ts), max(ts)))
plt.plot(trend, label='趋势')
# 生成网格
plt.grid()
plt.legend(loc='best',fontsize = 60)
plt.savefig(tqmc+'--趋势因素.png',bbox_inches='tight')
plt.subplot(413)
plt.plot(seasonal,label='季节')
# 生成网格
plt.grid()
plt.legend(loc='best',fontsize = 60)
plt.savefig(tqmc+'--季节因素.png',bbox_inches='tight')
plt.subplot(414)
plt.plot(residual, label='随机')
# 生成网格
plt.grid()
plt.legend(loc='best',fontsize = 60)
plt.savefig(tqmc+'--随机因素.png',bbox_inches='tight')
 plt.tight_layout()时间序列进行季节性分解之后的结果
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
rcParams['xtick.labelsize'] = 30
rcParams['ytick.labelsize'] = 30
#设置横纵坐标的名称以及对应字体格式
plt.tick_params(axis='both', which='major', labelsize=40)
plt.tick_params(axis='both', which='minor', labelsize=40)
def forecast_plt(tqmc,history,forecast):
plt.figure(figsize=(50,22))
plt.plot(history, label='历史负荷')
print(forecast.iloc[0])
plt.plot(forecast.iloc[0], color='red', label='预测曲线上界')
plt.plot(forecast.iloc[1], color='black', label='预测曲线均值',linestyle='dashed',linewidth=4)
plt.plot(forecast.iloc[2], color='blue', label='预测曲线下界',linestyle=':',linewidth=4)
plt.legend(loc='best',fontsize=30)
plt.title(tqmc+'历史负荷及预测结果',fontsize=30)
plt.savefig('原序列与预测序列--'+tqmc+'.png',bbox_inches='tight')
 plt.show()历史数据与总体预测结果
import datetime
#加上文件的生成日期，区别不同的文件版本
today = datetime.date.today().strftime('%m%d')
from statsmodels.tsa.seasonal import seasonal_decompose
for bs in industry_type_bs[0:15]:
tqmc = industry_type.loc[bs]['台区名称']
print(tqmc)
bs = str(bs)
try:
ts = data_time_series[bs]
except KeyError as e:
print(e)
continue
#如果负荷序列全为零值或者2018年的负荷序列都为零值，则跳过预测
if ~ts.any() or ~ts['2018'].any():
continue
ts = ts.apply(pd.to_numeric)
ts = ts.resample('W').max()
print(bs)
decomposition = seasonal_decompose(ts,extrapolate_trend='freq')
forecast_decomposition_plt(tqmc,ts,decomposition)
# 显式调用seasonal_decompose函数的extrapolate_trend参数可以强制提取残差，同时使趋势曲线更加平滑
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
residual = residual.fillna(method='ffill')
#如果时间序列的残差不平稳，则在预测之前要进行额外的一次差分操作
if ~test_stationarity(residual):
d = 1
else:
d = 0
#趋势预测
#y 是最近一年的负荷增长率
y = ts[str(ts.index.year[-1])].max() / ts[str(ts.index.year[-1]-1)].max() - 1
print(y)
#t是目前时间点位于的“S”型曲线的阶段
t = SCurveStage(y)
print(t)
trend_f = SCurveForecast(t,trend)
coefficent = proper_model(residual,20)
print(coefficent)
residual_f =compare_ARIMA_modes(residual,(coefficent[1], d, coefficent[2]))
if residual_f.shape[0]==0:
continue
seasonal_f = season_forcast_year(seasonal)
residual_f = residual_f[1:]
#三条预测曲线
forcast_model = trend_f+residual_f+seasonal_f
forcast_model_least = trend_f+seasonal_f
forcast_model_equal = (forcast_model+forcast_model_least)/2
df = pd.concat([forcast_model,forcast_model_equal,forcast_model_least],axis=1)
#画图
forecast_plt(tqmc,ts,df.T)
df.columns = [bs+'upper',bs,bs+'lower']
df = df.resample('M').max()
df = df.T
df_forecast = pd.concat([df_forecast,df])
#ruralcity 指城中村类别
 df_forecast.to_excel('Time_Series_Forecast_ruralcity_Monthly'+today+'.xlsx')只查看预测结果