更新时间：2022年4月21日

相关链接

（1）【第十届“泰迪杯”数据挖掘挑战赛】B题：电力系统负荷预测分析 问题一Baseline方案

（2）【第十届“泰迪杯”数据挖掘挑战赛】B题：电力系统负荷预测分析 问题一ARIMA、AutoARIMA、LSTM、Prophet 多方案实现

（3)【第十届“泰迪杯”数据挖掘挑战赛】B题：电力系统负荷预测分析 问题二 时间突变分析 Python实现

（4）【第十届“泰迪杯”数据挖掘挑战赛】B题：电力系统负荷预测分析 31页省一等奖论文及代码

完整代码下载链接

https://www.betterbench.top/#/35/detail

思路是设定一个阈值，大于或者小于这个阈值就判定为突变。阈值的设定有自定义或者使用统计方法进行设置。分析过程如下

1 定义绘图函数

import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
import plotly
import plotly.graph_objects as go
import plotly.io as pio
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from chart_studio import plotly as py
import matplotlib.pyplot as plt
import plotly.graph_objs as go
init_notebook_mode(connected=True)
import warnings
warnings.filterwarnings('ignore')
from pylab import mpl 
mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定默认字体 
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题

def plot_anomaly(ts,anomaly_pred = None,anomaly_true=None,file_name = 'file'):
    fig = go.Figure()
    yhat = go.Scatter(
    x = ts.index,
    y = ts,
    mode = 'lines',  name = ts.name)
    fig.add_trace(yhat)
    if  anomaly_pred is not None:
        status = go.Scatter(
        。。。略
    if  anomaly_true is not None:
        status = go.Scatter(
        x = anomaly_true.index,
        y = ts.loc[anomaly_true.index],
        。。。略
    fig.show()
    
def plot_anomaly_window(ts,anomaly_pred = None,file_name = 'file',window='1h'):
    fig = go.Figure()
    yhat = go.Scatter(
    x = ts.index,
    y = ts,
    。。。略
    fig.show()

2 通过对原始测量应用阈值来查找异常值

df = pd.read_csv('./data/附件1-区域15分钟负荷数据.csv')
df = df.rename(columns={'数据时间':'Time'})
df['Time']= pd.to_datetime(df['Time'])
df.set_index('Time',inplace=True)
col = df.columns[0]
df.plot()

df 

3 手动设置阈值

min_t = 203505
max_t = 300392.7345
df[col+'threshold_alarm'] = (df[col].clip(lower = min_t,upper=max_t) != df[col])
plot_anomaly(df[col],anomaly_pred = df[df[col+'threshold_alarm']==True][col+'threshold_alarm'],anomaly_true=None,file_name = 'file')

4 使用分位数设置阈值

min_t = df[col].quantile(0.03)
max_t = df[col].quantile(0.97)
df[col+'threshold_alarm'] = (df[col].clip(lower = min_t,upper=max_t) != df[col])
plot_anomaly(df[col],anomaly_pred = df[df[col+'threshold_alarm']==True][col+'threshold_alarm'],anomaly_true=None,file_name = 'file')

5 3Sigma原则（IQR）

[Q1−c×IQR，Q3+c×IQR]

Q1 = df[col].quantile(0.25)
Q3 = df[col].quantile(0.75)
IQR = Q3- Q1
。。。略


df[col+'threshold_alarm'] = (df[col].clip(lower = min_t,upper=max_t) != df[col])
plot_anomaly(df[col],anomaly_pred = df[df[col+'threshold_alarm']==True][col+'threshold_alarm'],anomaly_true=None,file_name = 'file')

6 设定分位数

window = 1
df[col+'_diff'] = df[col].diff(periods= window).fillna(0)

Q1 = df[col+'_diff'].quantile(0.25)
Q3 = df[col+'_diff'].quantile(0.75)
IQR = Q3- Q1
c = 2
min_t = Q1 - c*IQR
max_t = Q3 + c*IQR
df[col+'diff_alarm'] = (df[col+'_diff'].clip(lower = min_t,upper=max_t)!= df[col+'_diff'])
plot_anomaly_window(df[col],anomaly_pred = df[df[col+'diff_alarm']==True][col+'diff_alarm'],file_name = 'file',window=f'{window}h')

6.1 导入数据

import pandas as pd 
df = pd.read_csv('附件2-行业日负荷数据.csv')
metrics_df=pd.pivot_table(df,values='有功功率最大值（kw）',index='数据时间',columns='行业类型')
metrics_df.head()

metrics_df.reset_index(inplace=True)
metrics_df.fillna(0,inplace=True)
metrics_df.head()
metrics_df = metrics_df.rename(columns={'数据时间':'Time'})

df = metrics_df
df['Time']= pd.to_datetime(df['Time'])
df.set_index('Time',inplace=True)
col = df.columns[0]
df.plot()

col = df.columns[1]

Q1 = df[col].quantile(0.25)
Q3 = df[col].quantile(0.75)
IQR = Q3- Q1
。。。略

df[col+'threshold_alarm'] = (df[col].clip(lower = min_t,upper=max_t) != df[col])
plot_anomaly_window(df[col],anomaly_pred = df[df[col+'threshold_alarm']==True][col+'threshold_alarm'],window='1h',file_name = 'file')

7 根据分位数找到时间突变

window = 5
df[col+'ma'] = df[col].rolling(window=window,closed='left').mean()
kpi_col = col+'ma'+'diff'
df[kpi_col] = (df[col]-df[col+'ma']).fillna(0)


Q1 = df[kpi_col].quantile(0.25)
Q3 = df[kpi_col].quantile(0.75)
IQR = Q3- Q1
c = 2
min_t = Q1 - c*IQR
max_t = Q3 + c*IQR
df[kpi_col+'threshold_alarm'] = (df[kpi_col].clip(lower = min_t,upper=max_t) != df[kpi_col])
plot_anomaly_window(df[col],anomaly_pred = df[df[kpi_col+'threshold_alarm']==True][kpi_col+'threshold_alarm'],file_name = 'file',window=f'{window}h')