這是我之前工作做的一個項目

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')只查看預測結果

總結

以上是生活随笔為你收集整理的python用电度数设计_Python时间序列预测实战（电力负荷预测）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。