电力负荷案例之扩展思路

# 当前
feature_data['prev_1'] = feature_data['power_load'].shift(1)
feature_data['prev_2'] = feature_data['power_load'].shift(2)
feature_data['prev_3'] = feature_data['power_load'].shift(3)

# 优化后
for i in range(1, 25):
    feature_data[f'prev_{i}'] = feature_data['power_load'].shift(i)

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import logging
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error

# 解决中文乱码
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 生成模拟数据（含外部特征）
print("=== 1. 生成模拟数据（含外部特征） ===")

def generate_power_load_data_with_external(n_days=30):
    """生成电力负荷数据（含外部特征）"""
    np.random.seed(42)
    
    dates = pd.date_range('2024-01-01', periods=n_days*24, freq='h')
    power_load = []
    temperatures = []
    humidities = []
    wind_speeds = []
    
    for date in dates:
        hour = date.hour
        dayofweek = date.dayofweek
        month = date.month
        
        # 基础负荷
        if 0 <= hour < 6:
            base_load = 500
        elif 6 <= hour < 9:
            base_load = 800
        elif 9 <= hour < 17:
            base_load = 1000
        elif 17 <= hour < 21:
            base_load = 1200
        else:
            base_load = 700
        
        if dayofweek >= 5:
            base_load *= 0.8
        
        if month in [12, 1, 2]:
            base_load *= 1.2
        elif month in [6, 7, 8]:
            base_load *= 1.1
        
        # 模拟温度（季节性）
        if month in [12, 1, 2]:
            temp = np.random.normal(-5, 5)
        elif month in [6, 7, 8]:
            temp = np.random.normal(30, 5)
        else:
            temp = np.random.normal(15, 5)
        
        # 温度对负荷的影响
        if temp > 25:
            base_load *= 1.2  # 高温开空调
        elif temp < 0:
            base_load *= 1.1  # 低温开暖气
        
        # 模拟湿度
        humidity = np.random.uniform(30, 80)
        
        # 模拟风速
        wind_speed = np.random.uniform(0, 10)
        
        load = base_load + np.random.normal(0, 50)
        power_load.append(max(100, load))
        temperatures.append(temp)
        humidities.append(humidity)
        wind_speeds.append(wind_speed)
    
    data = pd.DataFrame({
        'time': dates,
        'power_load': power_load,
        'temperature': temperatures,
        'humidity': humidities,
        'wind_speed': wind_speeds
    })
    
    return data

data = generate_power_load_data_with_external(n_days=30)
print(f"数据形状: {data.shape}")
print(f"\n前5行数据:")
print(data.head())

# 2. 扩展特征工程
print("\n=== 2. 扩展特征工程 ===")

def extended_feature_engineering(data, window_size=24, logger=None):
    """
    扩展特征工程
    Args:
        data: 数据源
        window_size: 窗口大小
        logger: 日志对象
    Returns:
        feature_data: 特征数据
    """
    if logger:
        logger.info("开始扩展特征工程")
    
    feature_data = data.copy()
    
    # 1. 提取hour和month
    feature_data['hour'] = feature_data['time'].dt.hour
    feature_data['month'] = feature_data['time'].dt.month
    
    # 2. One-Hot编码
    feature_data = pd.get_dummies(feature_data, columns=['hour', 'month'])
    
    # 3. 添加窗口字段（扩展到24小时）
    for i in range(1, window_size + 1):
        feature_data[f'prev_{i}'] = feature_data['power_load'].shift(i)
    
    # 4. 获取昨日同时刻负荷
    feature_data['yesterday_time'] = feature_data['time'].apply(
        lambda x: (pd.to_datetime(x) - pd.Timedelta(days=1)).strftime('%Y-%m-%d %H:%M:%S')
    )
    
    time_load_map = dict(zip(
        feature_data['time'].dt.strftime('%Y-%m-%d %H:%M:%S'),
        feature_data['power_load']
    ))
    
    feature_data['yesterday_load'] = feature_data['yesterday_time'].map(time_load_map)
    
    # 5. 添加外部特征
    # 温度的滞后特征
    feature_data['temp_prev_1'] = feature_data['temperature'].shift(1)
    feature_data['temp_prev_2'] = feature_data['temperature'].shift(2)
    feature_data['temp_prev_3'] = feature_data['temperature'].shift(3)
    
    # 湿度的滞后特征
    feature_data['humidity_prev_1'] = feature_data['humidity'].shift(1)
    
    # 风速的滞后特征
    feature_data['wind_prev_1'] = feature_data['wind_speed'].shift(1)
    
    # 6. 添加统计特征
    # 过去24小时的统计量
    feature_data['load_mean_24'] = feature_data['power_load'].rolling(window=24).mean()
    feature_data['load_std_24'] = feature_data['power_load'].rolling(window=24).std()
    feature_data['load_max_24'] = feature_data['power_load'].rolling(window=24).max()
    feature_data['load_min_24'] = feature_data['power_load'].rolling(window=24).min()
    
    # 7. 删除NaN
    feature_data = feature_data.dropna().reset_index(drop=True)
    
    if logger:
        logger.info(f"扩展特征工程完成，形状: {feature_data.shape}")
    
    return feature_data

feature_data = extended_feature_engineering(data, window_size=24)
print(f"特征数据形状: {feature_data.shape}")
print(f"\n特征列名:")
exclude_columns = ['time', 'power_load', 'yesterday_time', 'temperature', 'humidity', 'wind_speed']
features = [col for col in feature_data.columns if col not in exclude_columns]
print(f"特征数量: {len(features)}")

# 3. 训练扩展模型
print("\n=== 3. 训练扩展模型 ===")

def train_extended_model(feature_data, features):
    """训练扩展模型"""
    # 数据集切分
    X = feature_data[features]
    y = feature_data['power_load']
    
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42, shuffle=False
    )
    
    # 训练模型
    model = GradientBoostingRegressor(
        n_estimators=100,
        max_depth=6,
        learning_rate=0.1,
        random_state=42
    )
    model.fit(X_train, y_train)
    
    # 评估
    y_pred = model.predict(X_test)
    mae = mean_absolute_error(y_test, y_pred)
    mse = mean_squared_error(y_test, y_pred)
    rmse = np.sqrt(mse)
    
    print(f"训练集: {X_train.shape}")
    print(f"测试集: {X_test.shape}")
    print(f"\n评估指标:")
    print(f"  MAE: {mae:.4f}")
    print(f"  MSE: {mse:.4f}")
    print(f"  RMSE: {rmse:.4f}")
    
    return model

model = train_extended_model(feature_data, features)

# 4. 特征重要性分析
print("\n=== 4. 特征重要性分析 ===")

def analyze_feature_importance(model, features, top_n=20):
    """分析特征重要性"""
    importances = model.feature_importances_
    
    # 排序
    indices = np.argsort(importances)[::-1]
    
    print(f"前{top_n}个重要特征:")
    for i, idx in enumerate(indices[:top_n]):
        print(f"  {i+1}. {features[idx]}: {importances[idx]:.4f}")
    
    # 可视化
    plt.figure(figsize=(12, 8))
    plt.bar(range(top_n), importances[indices[:top_n]], color='skyblue')
    plt.xlabel('特征')
    plt.ylabel('重要性')
    plt.title(f'前{top_n}个重要特征')
    plt.xticks(range(top_n), [features[i] for i in indices[:top_n]], rotation=90)
    plt.tight_layout()
    plt.savefig('extended_feature_importance.png', dpi=100)
    print(f"\n特征重要性图已保存")

analyze_feature_importance(model, features, top_n=20)

# 5. 算法对比
print("\n=== 5. 算法对比 ===")

def compare_algorithms(feature_data, features):
    """对比不同算法"""
    from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
    from sklearn.linear_model import LinearRegression
    
    # 数据集切分
    X = feature_data[features]
    y = feature_data['power_load']
    
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42, shuffle=False
    )
    
    algorithms = {
        'Linear Regression': LinearRegression(),
        'Random Forest': RandomForestRegressor(n_estimators=100, random_state=42),
        'Gradient Boosting': GradientBoostingRegressor(
            n_estimators=100, max_depth=6, learning_rate=0.1, random_state=42
        )
    }
    
    results = {}
    
    for name, model in algorithms.items():
        # 训练
        model.fit(X_train, y_train)
        
        # 预测
        y_pred = model.predict(X_test)
        
        # 评估
        mae = mean_absolute_error(y_test, y_pred)
        mse = mean_squared_error(y_test, y_pred)
        rmse = np.sqrt(mse)
        
        results[name] = {
            'MAE': mae,
            'MSE': mse,
            'RMSE': rmse
        }
        
        print(f"{name}:")
        print(f"  MAE: {mae:.4f}")
        print(f"  MSE: {mse:.4f}")
        print(f"  RMSE: {rmse:.4f}")
        print()
    
    return results

results = compare_algorithms(feature_data, features)

# 6. 扩展思路详解
print("\n=== 6. 扩展思路详解 ===")

def extension_ideas():
    """扩展思路详解"""
    print("""
扩展思路详解:

1. 特征工程优化
   1.1 时间窗口扩展
       - 现状: 最近3小时
       - 优化: 最近12小时、24小时
   
   1.2 历史负荷特征
       - 现状: 40个字段
       - 优化: 添加统计特征（均值、标准差、最大值、最小值）
   
   1.3 外部特征
       - 温度
       - 湿度
       - 风速
       - 历史工业用电量
       - 历史居民用电量

2. 算法优化
   2.1 尝试不同算法
       - XGBoost
       - LightGBM
       - CatBoost
       - 神经网络（LSTM、GRU）
   
   2.2 集成学习
       - 多模型集成
       - 加权平均
       - Stacking

3. 预测优化
   3.1 滚动预测
       - 每次预测后更新数据
       - 滚动窗口
   
   3.2 多步预测
       - 预测未来多个时间点
       - 直接多步预测
       - 递归多步预测
   
   3.3 集成预测
       - 多模型集成预测
       - 提高稳定性

4. 数据来源
   4.1 天聚地合
       - 网址: 天聚地合
       - 数据: 股票、天气、空气质量
       - 免费: 每天5次
       - 付费: 不限次数
   
   4.2 其他数据源
       - 天气预报
       - 历史天气
       - 空气质量
""")

extension_ideas()

# 7. 外部特征详解
print("\n=== 7. 外部特征详解 ===")

def external_features_explanation():
    """外部特征详解"""
    print("""
外部特征详解:

1. 温度
   - 影响: 高温开空调，低温开暖气
   - 获取: 天气预报、历史天气
   - 特征: 当前温度、滞后温度

2. 湿度
   - 影响: 影响体感温度
   - 获取: 天气预报
   - 特征: 当前湿度、滞后湿度

3. 风速
   - 影响: 影响风力发电
   - 获取: 天气预报
   - 特征: 当前风速、滞后风速

4. 历史用电量
   - 工业: 工业用电量
   - 居民: 居民用电量
   - 商业: 商业用电量

5. 节假日
   - 影响: 节假日用电模式不同
   - 获取: 节假日API
   - 特征: 是否节假日、节假日类型
""")

external_features_explanation()

# 8. 数据平台介绍
print("\n=== 8. 数据平台介绍 ===")

def data_platform_introduction():
    """数据平台介绍"""
    print("""
数据平台介绍:

1. 天聚地合
   1.1 功能
       - 提供各种数据集
       - 定制化数据服务
       - API接口
   
   1.2 数据类型
       - 股票数据
       - 天气数据
       - 空气质量数据
       - 银行数据
       - 身份证数据
       - 基站数据
   
   1.3 收费
       - 免费: 每天5次
       - 付费: 黑钻最高不限次数
   
   1.4 示例
       - 纽约股票交易所数据
       - 标准普尔500数据
       - 上证指数数据
       - 全国城市空气质量
       - 历史天气数据

2. 其他数据源
   2.1 天气预报
       - 越来越精准
       - 可以获取未来天气
   
   2.2 历史天气
       - 用于训练模型
       - 分析历史趋势
   
   2.3 实时数据
       - 爬虫获取
       - 注意合法性
""")

data_platform_introduction()

# 9. 完整扩展代码
print("\n=== 9. 完整扩展代码 ===")

def complete_extension_code():
    """完整扩展代码"""
    print("""
# 1. 扩展特征工程
def extended_feature_engineering(data, window_size=24):
    feature_data = data.copy()
    
    # 提取hour和month
    feature_data['hour'] = feature_data['time'].dt.hour
    feature_data['month'] = feature_data['time'].dt.month
    
    # One-Hot编码
    feature_data = pd.get_dummies(feature_data, columns=['hour', 'month'])
    
    # 添加窗口字段（扩展到24小时）
    for i in range(1, window_size + 1):
        feature_data[f'prev_{i}'] = feature_data['power_load'].shift(i)
    
    # 获取昨日同时刻负荷
    feature_data['yesterday_time'] = feature_data['time'].apply(
        lambda x: (pd.to_datetime(x) - pd.Timedelta(days=1)).strftime('%Y-%m-%d %H:%M:%S')
    )
    time_load_map = dict(zip(
        feature_data['time'].dt.strftime('%Y-%m-%d %H:%M:%S'),
        feature_data['power_load']
    ))
    feature_data['yesterday_load'] = feature_data['yesterday_time'].map(time_load_map)
    
    # 添加外部特征
    feature_data['temp_prev_1'] = feature_data['temperature'].shift(1)
    feature_data['humidity_prev_1'] = feature_data['humidity'].shift(1)
    feature_data['wind_prev_1'] = feature_data['wind_speed'].shift(1)
    
    # 添加统计特征
    feature_data['load_mean_24'] = feature_data['power_load'].rolling(window=24).mean()
    feature_data['load_std_24'] = feature_data['power_load'].rolling(window=24).std()
    feature_data['load_max_24'] = feature_data['power_load'].rolling(window=24).max()
    feature_data['load_min_24'] = feature_data['power_load'].rolling(window=24).min()
    
    # 删除NaN
    feature_data = feature_data.dropna().reset_index(drop=True)
    
    return feature_data

# 2. 算法对比
def compare_algorithms(X_train, X_test, y_train, y_test):
    algorithms = {
        'Linear Regression': LinearRegression(),
        'Random Forest': RandomForestRegressor(n_estimators=100),
        'Gradient Boosting': GradientBoostingRegressor(n_estimators=100),
        'XGBoost': XGBRegressor(n_estimators=100)
    }
    
    for name, model in algorithms.items():
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        mae = mean_absolute_error(y_test, y_pred)
        print(f"{name}: MAE = {mae:.4f}")
""")

complete_extension_code()

# 10. 总结
def extension_summary():
    """扩展思路总结"""
    print("=" * 50)
    print("扩展思路总结")
    print("=" * 50)
    
    print("\n1. 社会价值")
    print("   预测误差越小，节省的钱越多")
    
    print("\n2. 特征工程优化")
    print("   时间窗口扩展")
    print("   历史负荷特征")
    print("   外部特征（温度、湿度、风速）")
    
    print("\n3. 算法优化")
    print("   尝试不同算法")
    print("   集成学习")
    
    print("\n4. 预测优化")
    print("   滚动预测")
    print("   多步预测")
    print("   集成预测")
    
    print("\n5. 数据来源")
    print("   天聚地合")
    print("   天气预报")
    print("   历史天气")
    
    print("\n" + "=" * 50)
    print("扩展思路完成!")
    print("=" * 50)

extension_summary()