import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.data import TensorDataset, DataLoader
from sklearn.preprocessing import MinMaxScaler
from vmdpy import VMD
from sklearn.cluster import OPTICS
import math
import matplotlib.pyplot as plt
from torch.cuda import amp
import time
import logging
# 设置日志
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger(__name__)
# -------------------- 配置参数 --------------------
CONFIG = {
'seq_length': 12, # 输入序列长度
'batch_size': 256, # 批大小
'epochs': 150, # 训练轮数
'd_model': 64, # 模型维度
'nhead': 4, # 注意力头数
'num_layers': 2, # Transformer层数
'dim_feedforward': 128, # 前馈网络维度
'dropout': 0.1, # Dropout率
'vmd_k': 8, # VMD分解模态数
'vmd_alpha': 2000, # VMD带宽约束
'vmd_tau': 0.1, # VMD噪声容忍度
'sru_hidden': 32, # SRU隐藏单元数
'sru_layers': 2, # SRU层数
'cluster_min_samples': 3, # 聚类最小样本数
'use_amp': True, # 启用混合精度训练
'device': 'cuda' if torch.cuda.is_available() else 'cpu'
}
# -------------------- 数据预处理模块 --------------------
def load_data(file_path):
"""加载并预处理风电数据"""
logger.info(f"加载数据: {file_path}")
try:
data = pd.read_csv(file_path)
logger.info(f"数据列: {data.columns.tolist()}")
except Exception as e:
logger.error(f"加载文件错误: {e}")
raise
time_col = 'Time'
power_col = 'Power'
# 检查必要列是否存在
required_cols = [time_col, power_col]
missing_cols = [col for col in required_cols if col not in data.columns]
if missing_cols:
logger.warning(f"警告: 缺少必要列 {missing_cols}")
# 如果缺少Time列，创建一个简单的时间索引
if time_col not in data.columns:
logger.info("创建简单时间索引")
data[time_col] = pd.date_range(start='2023-01-01', periods=len(data), freq='H')
# 时间序列处理
data[time_col] = pd.to_datetime(data[time_col])
data = data.set_index(time_col).resample('1h').asfreq()
data.index = data.index.fillna(pd.Timestamp.now())
# 功率处理
if power_col not in data.columns:
logger.error("错误: 缺少功率列")
raise ValueError("CSV文件必须包含功率列")
data[power_col] = data[power_col].replace(0, np.nan)
data[power_col] = data[power_col].interpolate(method='time')
data[power_col] = data[power_col].fillna(0)
# 添加风速特征 - 使用 windspeed_10m 和 windspeed_100m
if 'windspeed_10m' in data.columns and 'windspeed_100m' in data.columns:
logger.info("使用10m和100m风速数据计算有效风速")
# 对于NaN值，进行插值
data['windspeed_10m'] = data['windspeed_10m'].interpolate(method='time').fillna(method='bfill')
data['windspeed_100m'] = data['windspeed_100m'].interpolate(method='time').fillna(method='bfill')
# 计算平均风速作为有效风速
data['EffectiveWind'] = (data['windspeed_10m'] + data['windspeed_100m']) / 2
else:
logger.warning("风速数据缺失，生成模拟风速数据")
# 基于功率生成合理的模拟风速
max_power = data[power_col].max()
# 基本风速曲线：功率越高风速越高（但非线性）
base_wind = np.clip(data[power_col] / max_power * 10, 2, 25)
# 添加随机波动
random_fluctuation = np.random.normal(0, 1.5, len(data))
data['EffectiveWind'] = base_wind + random_fluctuation
# 确保只返回必要的列
result_cols = ['Power']
if 'EffectiveWind' in data.columns:
result_cols.append('EffectiveWind')
logger.info(f"返回列: {result_cols}")
return data[result_cols].values
# -------------------- VMD分解模块（替代CEEMDAN）-------------------
def vmd_decomposition(data, alpha=2000, tau=0.1, K=8):
"""执行VMD分解 - 比CEEMDAN更快更高效"""
logger.info("开始VMD分解...")
start_time = time.time()
# 仅对功率数据进行分解
power_data = data[:, 0].flatten()
# 执行VMD分解
u, u_hat, omega = VMD(power_data, alpha, tau, K, DC=0, init=1, tol=1e-7)
# 保留分解结果和气象特征
imfs = u
residual = power_data - np.sum(imfs, axis=0)
# 添加残差作为最后一个分量
components = np.vstack([imfs, residual])
# 将气象特征附加到每个分量
components_with_features = []
for i in range(len(components)):
comp = components[i]
# 为每个分量添加气象特征
# 如果数据只有一列（功率），则复制该列作为特征
if data.shape[1] == 1:
comp_with_features = np.column_stack((comp, np.zeros_like(comp)))
else:
comp_with_features = np.column_stack((comp, data[:, 1]))
components_with_features.append(comp_with_features)
logger.info(f"VMD分解完成，耗时: {time.time() - start_time:.2f}秒")
return components_with_features
# -------------------- 多维度特征提取 --------------------
def extract_features(component):
"""提取序列的多维度特征"""
power_series = component[:, 0]
# 1. 样本熵
def sample_entropy(series, m=2, alpha=0.2):
n = len(series)
if n < m + 1: return 0
std = np.std(series)
r = alpha * std
def _phi(_m):
x = np.array([series[i:i + _m] for i in range(n - _m + 1)])
C = 0
for i in range(len(x)):
dist = np.max(np.abs(x[i] - x), axis=1)
C += np.sum((dist < r) & (dist > 0))
return C / ((n - _m) * (n - _m + 1))
return -np.log(_phi(m + 1) / _phi(m)) if _phi(m) != 0 else 0
# 2. 排列熵
def permutation_entropy(series, d=3, tau=1):
n = len(series)
if n < d * tau: return 0
# 创建符号序列
permutations = []
for i in range(n - d * tau + 1):
segment = series[i:i + d * tau:tau]
permutations.append(tuple(np.argsort(segment)))
# 计算概率分布
unique, counts = np.unique(permutations, return_counts=True)
probs = counts / len(permutations)
# 计算熵
return -np.sum(probs * np.log(probs))
# 3. 频域能量
fft_vals = np.abs(np.fft.rfft(power_series))
spectral_energy = np.sum(fft_vals[:len(fft_vals) // 2]) / np.sum(fft_vals)
return np.array([
sample_entropy(power_series),
permutation_entropy(power_series),
spectral_energy
])
# -------------------- 轻量化Transformer模型（高频序列）-------------------
class ProbSparseAttention(nn.Module):
"""概率稀疏注意力机制 - 降低计算复杂度"""
def __init__(self, d_model, n_heads, factor=5):
super().__init__()
self.d_model = d_model
self.n_heads = n_heads
self.factor = factor
self.head_dim = d_model // n_heads
def forward(self, Q, K, V):
batch_size, seq_len, _ = Q.size()
# 采样关键点
M = self.factor * int(np.ceil(np.log(seq_len)))
sample_indices = torch.randperm(seq_len)[:M]
K_sampled = K[:, sample_indices, :]
V_sampled = V[:, sample_indices, :]
# 计算稀疏注意力
Q = Q.view(batch_size, seq_len, self.n_heads, self.head_dim).transpose(1, 2)
K_sampled = K_sampled.view(batch_size, M, self.n_heads, self.head_dim).transpose(1, 2)
V_sampled = V_sampled.view(batch_size, M, self.n_heads, self.head_dim).transpose(1, 2)
attn_scores = torch.matmul(Q, K_sampled.transpose(-2, -1)) / np.sqrt(self.head_dim)
attn_weights = F.softmax(attn_scores, dim=-1)
output = torch.matmul(attn_weights, V_sampled)
output = output.transpose(1, 2).contiguous().view(batch_size, seq_len, self.d_model)
return output
class PositionalEncoding(nn.Module):
def __init__(self, d_model, dropout=0.1, max_len=5000):
super().__init__()
self.dropout = nn.Dropout(p=dropout)
pe = torch.zeros(max_len, d_model)
position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1)
div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-math.log(10000.0) / d_model))
pe[:, 0::2] = torch.sin(position * div_term)
pe[:, 1::2] = torch.cos(position * div_term)
pe = pe.unsqueeze(0)
self.register_buffer('pe', pe)
def forward(self, x):
x = x + self.pe[:, :x.size(1)]
return self.dropout(x)
class DistillingLayer(nn.Module):
"""蒸馏层 - 压缩序列长度"""
def __init__(self, d_model):
super().__init__()
self.conv = nn.Conv1d(
in_channels=d_model,
out_channels=d_model,
kernel_size=3,
stride=2,
padding=1
)
self.activation = nn.ReLU()
def forward(self, x):
# x: [batch, seq_len, d_model]
x = x.permute(0, 2, 1) # [batch, d_model, seq_len]
x = self.conv(x)
x = self.activation(x)
return x.permute(0, 2, 1) # [batch, new_seq, d_model]
class EfficientTransformer(nn.Module):
"""高效Transformer模型 - 使用概率稀疏注意力和蒸馏机制"""
def __init__(self, input_dim, d_model=64, nhead=4, num_layers=2, dim_feedforward=128, dropout=0.1):
super().__init__()
self.d_model = d_model
self.embedding = nn.Linear(input_dim, d_model)
self.pos_encoder = PositionalEncoding(d_model, dropout)
# 编码器层
self.encoder_layers = nn.ModuleList()
for i in range(num_layers):
self.encoder_layers.append(nn.TransformerEncoderLayer(
d_model=d_model, nhead=nhead,
dim_feedforward=dim_feedforward,
dropout=dropout, batch_first=True
))
# 蒸馏层
self.distill_layers = nn.ModuleList([
DistillingLayer(d_model) for _ in range(num_layers - 1)
])
# 解码器
self.decoder = nn.TransformerDecoder(
nn.TransformerDecoderLayer(
d_model=d_model, nhead=nhead,
dim_feedforward=dim_feedforward, dropout=dropout,
batch_first=True
), num_layers=1
)
self.output_layer = nn.Linear(d_model, 1)
def forward(self, src, tgt=None):
# 嵌入和位置编码
src = self.embedding(src) * math.sqrt(self.d_model)
src = self.pos_encoder(src)
# 编码过程
for i, layer in enumerate(self.encoder_layers):
src = layer(src)
if i < len(self.distill_layers):
src = self.distill_layers[i](src)
# 解码过程
if tgt is None:
tgt = torch.zeros(src.size(0), src.size(1), self.d_model, device=src.device)
else:
tgt = self.embedding(tgt)
tgt = self.pos_encoder(tgt)
output = self.decoder(tgt, src)
output = self.output_layer(output)
return output[:, -1, :].squeeze(-1)
# -------------------- SRU模型（低频序列）-------------------
class SRU(nn.Module):
"""Simple Recurrent Unit - 比GRU更快的替代方案"""
def __init__(self, input_size, hidden_size, num_layers=2):
super().__init__()
self.hidden_size = hidden_size
self.num_layers = num_layers
# 门控参数
self.gates = nn.ModuleList()
for i in range(num_layers):
in_dim = input_size if i == 0 else hidden_size
self.gates.append(nn.Linear(in_dim + hidden_size, 3 * hidden_size))
def forward(self, x):
batch_size, seq_len, _ = x.size()
h = torch.zeros(self.num_layers, batch_size, self.hidden_size, device=x.device)
outputs = []
for t in range(seq_len):
input_t = x[:, t, :]
new_h = []
for i in range(self.num_layers):
# 当前层的输入
layer_input = input_t if i == 0 else new_h[i - 1]
# 与前一隐藏状态连接
combined = torch.cat((layer_input, h[i]), dim=1)
# 计算门控
gates = self.gates[i](combined)
f, r, c = torch.chunk(gates, 3, dim=1)
f = torch.sigmoid(f)
r = torch.sigmoid(r)
c = torch.tanh(c)
# 更新隐藏状态
h_i = f * h[i] + (1 - f) * c
output = r * h_i
new_h.append(h_i)
input_t = output
outputs.append(output)
h = torch.stack(new_h, dim=0)
return torch.stack(outputs, dim=1)
class SRUAttention(nn.Module):
"""带注意力机制的SRU模型"""
def __init__(self, input_dim, hidden_size=32, num_layers=2):
super().__init__()
self.sru = SRU(input_dim, hidden_size, num_layers)
self.attention = nn.Sequential(
nn.Linear(hidden_size, hidden_size),
nn.Tanh(),
nn.Linear(hidden_size, 1),
nn.Softmax(dim=1)
)
self.fc = nn.Linear(hidden_size, 1)
def forward(self, x):
# SRU输出: [batch, seq_len, hidden_size]
sru_out = self.sru(x)
# 注意力权重
attn_weights = self.attention(sru_out)
# 上下文向量
context = torch.sum(attn_weights * sru_out, dim=1)
return self.fc(context).squeeze(-1)
# -------------------- 模型路由网络 --------------------
class RoutingNetwork(nn.Module):
"""动态模型路由网络 - 根据序列特征选择模型"""
def __init__(self, transformer, sru_model):
super().__init__()
self.transformer = transformer
self.sru_model = sru_model
self.router = nn.Sequential(
nn.Linear(3, 16), # 输入特征数
nn.ReLU(),
nn.Linear(16, 2),
nn.Softmax(dim=1)
)
def forward(self, x, features):
# 特征: [样本熵, 排列熵, 频域能量]
route_probs = self.router(features)
# 使用两个模型进行预测
trans_pred = self.transformer(x)
sru_pred = self.sru_model(x)
# 加权组合
return (route_probs[:, 0] * trans_pred +
route_probs[:, 1] * sru_pred)
# -------------------- 主流程（含动态路由策略）-------------------
if __name__ == "__main__":
try:
# 设置随机种子确保可重复性
torch.manual_seed(42)
np.random.seed(42)
# 数据加载与预处理
file_path = 'G:/shuju/Location1.csv'
raw_data = load_data(file_path)
# 打印前5行数据
logger.info(f"数据形状: {raw_data.shape}")
logger.info(f"前5行数据:\n{raw_data[:5]}")
# 数据标准化
scalers = []
scaled_data = np.zeros_like(raw_data)
for i in range(raw_data.shape[1]):
scaler = MinMaxScaler(feature_range=(0, 1))
scaled_data[:, i] = scaler.fit_transform(raw_data[:, i].reshape(-1, 1)).flatten()
scalers.append(scaler)
# VMD分解（替代CEEMDAN）
components = vmd_decomposition(scaled_data,
alpha=CONFIG['vmd_alpha'],
tau=CONFIG['vmd_tau'],
K=CONFIG['vmd_k'])
# 提取特征并聚类分组
logger.info("提取特征并聚类分组...")
features = np.array([extract_features(comp) for comp in components])
# OPTICS聚类
clusterer = OPTICS(min_samples=CONFIG['cluster_min_samples'], xi=0.05)
labels = clusterer.fit_predict(features)
# 按聚类结果分组
grouped_components = {}
for i, label in enumerate(labels):
if label not in grouped_components:
grouped_components[label] = []
grouped_components[label].append(components[i][:, 0]) # 只取功率部分
# 重构序列
reconstructed_series = []
for label, comp_list in grouped_components.items():
if len(comp_list) > 1:
reconstructed_series.append(np.sum(comp_list, axis=0))
else:
reconstructed_series.append(comp_list[0])
logger.info(f"重构为 {len(reconstructed_series)} 个序列")
# 转换为监督学习格式
def create_sequences(data, seq_length):
X, y = [], []
for i in range(len(data) - seq_length):
X.append(data[i:i + seq_length])
y.append(data[i + seq_length])
return np.array(X), np.array(y)
# 创建数据集
datasets = []
for series in reconstructed_series:
X, y = create_sequences(series, CONFIG['seq_length'])
datasets.append((X, y))
# 模型初始化
transformer = EfficientTransformer(input_dim=1,
d_model=CONFIG['d_model'],
nhead=CONFIG['nhead'],
num_layers=CONFIG['num_layers'],
dim_feedforward=CONFIG['dim_feedforward'],
dropout=CONFIG['dropout']).to(CONFIG['device'])
sru_model = SRUAttention(input_dim=1,
hidden_size=CONFIG['sru_hidden'],
num_layers=CONFIG['sru_layers']).to(CONFIG['device'])
routing_net = RoutingNetwork(transformer, sru_model).to(CONFIG['device'])
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(routing_net.parameters(), lr=0.001)
# 混合精度训练的梯度缩放器
scaler = amp.GradScaler(enabled=CONFIG['use_amp'])
# 训练循环
logger.info("开始训练...")
start_time = time.time()
# 由于我们有多个序列，需要合并训练数据
all_X, all_y, all_features = [], [], []
for i, (X, y) in enumerate(datasets):
# 为每个序列提取特征
seq_features = features[i]
all_X.append(X)
all_y.append(y)
# 为每个样本复制特征
all_features.extend([seq_features] * len(X))
all_X = np.concatenate(all_X)
all_y = np.concatenate(all_y)
all_features = np.array(all_features)
# 按时间顺序划分数据集
split_index = int(len(all_X) * 0.9)
train_dataset = TensorDataset(
torch.tensor(all_X[:split_index], dtype=torch.float32),
torch.tensor(all_y[:split_index], dtype=torch.float32),
torch.tensor(all_features[:split_index], dtype=torch.float32)
)
test_dataset = TensorDataset(
torch.tensor(all_X[split_index:], dtype=torch.float32),
torch.tensor(all_y[split_index:], dtype=torch.float32),
torch.tensor(all_features[split_index:], dtype=torch.float32)
)
# 创建DataLoader
train_loader = DataLoader(train_dataset, batch_size=CONFIG['batch_size'], shuffle=True)
test_loader = DataLoader(test_dataset, batch_size=CONFIG['batch_size'])
# 训练循环
for epoch in range(CONFIG['epochs']):
routing_net.train()
epoch_loss = 0
for inputs, targets, feat in train_loader:
inputs, targets, feat = inputs.to(CONFIG['device']), targets.to(CONFIG['device']), feat.to(
CONFIG['device'])
optimizer.zero_grad()
# 混合精度训练
with amp.autocast(enabled=CONFIG['use_amp']):
outputs = routing_net(inputs.unsqueeze(-1), feat)
loss = criterion(outputs, targets)
scaler.scale(loss).backward()
scaler.step(optimizer)
scaler.update()
epoch_loss += loss.item()
avg_loss = epoch_loss / len(train_loader)
logger.info(f"Epoch {epoch + 1}/{CONFIG['epochs']}, Loss: {avg_loss:.6f}")
logger.info(f"训练完成，总耗时: {time.time() - start_time:.2f}秒")
# 预测与评估
def inverse_transform(predictions, feature_idx=0):
return scalers[feature_idx].inverse_transform(predictions.reshape(-1, 1))
routing_net.eval()
all_preds, all_targets = [], []
with torch.no_grad():
for inputs, targets, feat in test_loader:
inputs, targets, feat = inputs.to(CONFIG['device']), targets.to(CONFIG['device']), feat.to(
CONFIG['device'])
outputs = routing_net(inputs.unsqueeze(-1), feat)
all_preds.append(outputs.cpu().numpy())
all_targets.append(targets.cpu().numpy())
all_preds = np.concatenate(all_preds)
all_targets = np.concatenate(all_targets)
# 反归一化
final_pred = inverse_transform(all_preds)
y_true = inverse_transform(all_targets)
# 评估指标
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
rmse = np.sqrt(mean_squared_error(y_true, final_pred))
mae = mean_absolute_error(y_true, final_pred)
r2 = r2_score(y_true, final_pred)
logger.info(f"最终评估 - RMSE: {rmse:.4f}, MAE: {mae:.4f}, R²: {r2:.4f}")
# 绘制预测结果与真实值的对比图
plt.figure(figsize=(15, 6))
plt.plot(y_true[:500], label='True', linewidth=2)
plt.plot(final_pred[:500], label='Predicted', linestyle='--')
plt.title(f'Wind Power Prediction\nRMSE: {rmse:.2f}, MAE: {mae:.2f}, R²: {r2:.4f}')
plt.xlabel('Time Steps')
plt.ylabel('Power (MW)')
plt.legend()
plt.grid(True)
plt.savefig('optimized_prediction_comparison.png', dpi=300)
plt.show()
# 保存模型
torch.save(routing_net.state_dict(), 'optimized_wind_power_model.pth')
logger.info("模型已保存")
except Exception as e:
logger.error(f"程序出错: {str(e)}")
import traceback
logger.error(traceback.format_exc())阅读该代码生成流程图，并确保流程图美观，简洁符合代码
本文介绍了一种使用递归方法实现的排列算法，通过交换数组中元素的位置来生成所有可能的排列组合，并提供了两种不同的C++代码实现方案。
扫一扫
79
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
