WGAN-GP与Relief算法在电力系统动态安全域边界生成中的应用

在电力系统动态安全分析中，快速生成准确的动态安全域边界对实时安全评估至关重要。传统逐点校验法虽然精度高，但计算成本巨大，难以满足在线应用需求。而基于WGAN-GP的区域法能够从有限样本中学习安全域边界的分布特征，实现边界的快速生成和动态更新。

本文面向电力系统分析工程师和机器学习研究者，将详细讲解如何结合WGAN-GP和Relief算法构建动态安全域边界生成模型。通过完整的代码实现和参数分析，展示该方法相比传统最小二乘法拟合的显著优势。

1. 理解动态安全域边界生成的核心问题

1.1 动态安全域的定义与重要性

动态安全域描述了电力系统在给定运行状态下，能够保持暂态稳定的关键参数空间边界。在实际系统中，这个边界通常是高维空间中的复杂曲面，传统方法需要大量时域仿真来验证边界点的安全性。

1.2 逐点校验法的局限性

逐点校验通过扫描参数空间并执行时域仿真来判定安全性，每个边界点都需要独立的仿真计算。对于n维参数空间，计算复杂度随维度指数增长，在实时安全评估中几乎不可行。

1.3 区域法的基本思想

区域法通过机器学习模型从有限样本中学习安全边界特征，然后快速预测新样本的安全性。WGAN-GP作为生成对抗网络的改进版本，能够稳定地学习复杂分布，适合安全域边界建模。

2. 环境准备与核心算法原理

2.1 环境依赖配置

项目需要Python 3.8+环境，主要依赖包包括：

# requirements.txt
torch>=1.9.0
numpy>=1.21.0
scikit-learn>=1.0.0
matplotlib>=3.5.0
pandas>=1.3.0
scipy>=1.7.0

关键依赖说明：

• PyTorch：WGAN-GP模型实现的基础框架
• scikit-learn：用于数据预处理和Relief特征选择
• scipy：提供最小二乘法拟合的对比实现

2.2 WGAN-GP算法原理

WGAN-GP（Wasserstein GAN with Gradient Penalty）通过Wasserstein距离衡量真实分布与生成分布的差异，解决了原始GAN训练不稳定的问题。

关键改进点：

• 使用Wasserstein距离替代JS散度，提供更有意义的训练信号
• 通过梯度惩罚项强制满足Lipschitz约束，避免权重裁剪带来的优化问题

损失函数定义：

def wgan_gp_loss(real_data, fake_data, discriminator, lambda_gp=10):
    # 计算基本Wasserstein损失
    real_loss = torch.mean(discriminator(real_data))
    fake_loss = torch.mean(discriminator(fake_data))
    w_loss = fake_loss - real_loss
    
    # 梯度惩罚项
    alpha = torch.rand(real_data.size(0), 1)
    interpolates = alpha * real_data + (1 - alpha) * fake_data
    interpolates.requires_grad_(True)
    d_interpolates = discriminator(interpolates)
    
    gradients = torch.autograd.grad(
        outputs=d_interpolates, inputs=interpolates,
        grad_outputs=torch.ones_like(d_interpolates),
        create_graph=True, retain_graph=True
    )[0]
    
    gradient_penalty = ((gradients.norm(2, dim=1) - 1) ** 2).mean()
    
    return w_loss + lambda_gp * gradient_penalty

2.3 Relief特征选择算法

Relief算法通过评估特征在同类和异类样本中的区分能力进行特征选择，特别适合高维电力系统参数的选择。

from sklearn.neighbors import NearestNeighbors

class ReliefSelector:
    def __init__(self, n_features_to_select=10):
        self.n_features = n_features_to_select
        
    def fit(self, X, y):
        n_samples, n_features = X.shape
        weights = np.zeros(n_features)
        
        # 找到每个样本的最近邻
        nn = NearestNeighbors(n_neighbors=2)
        nn.fit(X)
        
        for i in range(n_samples):
            # 找到同类的最近邻和异类的最近邻
            distances, indices = nn.kneighbors(X[i:i+1])
            hit_idx = indices[0][1]  # 最近邻（排除自身）
            
            # 找到异类最近邻
            other_class_mask = y != y[i]
            if np.any(other_class_mask):
                miss_idx = np.argmin(np.linalg.norm(
                    X[other_class_mask] - X[i], axis=1))
                miss_sample = X[other_class_mask][miss_idx]
                
                # 更新特征权重
                weights -= np.abs(X[i] - X[hit_idx]) / n_samples
                weights += np.abs(X[i] - miss_sample) / n_samples
                
        self.feature_weights = weights
        self.selected_features = np.argsort(weights)[-self.n_features:]
        
    def transform(self, X):
        return X[:, self.selected_features]

3. 完整实现：基于WGAN-GP的动态安全域边界生成

3.1 数据准备与预处理

电力系统动态安全数据通常包含发电机功角、电压、功率等参数，以及对应的安全标签（0/1表示不安全/安全）。

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

class SecurityDomainData:
    def __init__(self, data_path):
        self.data = pd.read_csv(data_path)
        self.scaler = StandardScaler()
        
    def preprocess(self, test_size=0.2):
        # 分离特征和标签
        X = self.data.drop('security_label', axis=1).values
        y = self.data['security_label'].values
        
        # 特征选择
        selector = ReliefSelector(n_features_to_select=15)
        X_selected = selector.fit_transform(X, y)
        
        # 数据标准化
        X_scaled = self.scaler.fit_transform(X_selected)
        
        # 划分训练测试集
        X_train, X_test, y_train, y_test = train_test_split(
            X_scaled, y, test_size=test_size, random_state=42)
        
        return X_train, X_test, y_train, y_test, selector.selected_features

3.2 WGAN-GP模型实现

import torch
import torch.nn as nn

class Generator(nn.Module):
    def __init__(self, latent_dim, output_dim):
        super(Generator, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(latent_dim, 128),
            nn.ReLU(),
            nn.Linear(128, 256),
            nn.ReLU(),
            nn.Linear(256, 512),
            nn.ReLU(),
            nn.Linear(512, output_dim),
            nn.Tanh()
        )
        
    def forward(self, z):
        return self.net(z)

class Discriminator(nn.Module):
    def __init__(self, input_dim):
        super(Discriminator, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, 512),
            nn.LeakyReLU(0.2),
            nn.Linear(512, 256),
            nn.LeakyReLU(0.2),
            nn.Linear(256, 128),
            nn.LeakyReLU(0.2),
            nn.Linear(128, 1)
        )
        
    def forward(self, x):
        return self.net(x)

class WGAN_GP:
    def __init__(self, latent_dim, data_dim, device='cuda'):
        self.generator = Generator(latent_dim, data_dim).to(device)
        self.discriminator = Discriminator(data_dim).to(device)
        self.device = device
        
        # 优化器设置
        self.g_optimizer = torch.optim.Adam(
            self.generator.parameters(), lr=0.0001, betas=(0.5, 0.9))
        self.d_optimizer = torch.optim.Adam(
            self.discriminator.parameters(), lr=0.0001, betas=(0.5, 0.9))
        
    def train_epoch(self, real_data, batch_size=64, n_critic=5):
        real_data = torch.FloatTensor(real_data).to(self.device)
        
        for _ in range(n_critic):
            # 训练判别器
            self.d_optimizer.zero_grad()
            
            # 生成假数据
            z = torch.randn(batch_size, self.generator.net[0].in_features).to(self.device)
            fake_data = self.generator(z)
            
            # 计算损失
            d_loss = self.compute_gradient_penalty(real_data, fake_data)
            d_loss.backward()
            self.d_optimizer.step()
        
        # 训练生成器
        self.g_optimizer.zero_grad()
        z = torch.randn(batch_size, self.generator.net[0].in_features).to(self.device)
        fake_data = self.generator(z)
        g_loss = -torch.mean(self.discriminator(fake_data))
        g_loss.backward()
        self.g_optimizer.step()
        
        return d_loss.item(), g_loss.item()
    
    def compute_gradient_penalty(self, real_data, fake_data):
        # WGAN-GP梯度惩罚计算
        alpha = torch.rand(real_data.size(0), 1).to(self.device)
        interpolates = (alpha * real_data + (1 - alpha) * fake_data).requires_grad_(True)
        d_interpolates = self.discriminator(interpolates)
        
        gradients = torch.autograd.grad(
            outputs=d_interpolates, inputs=interpolates,
            grad_outputs=torch.ones_like(d_interpolates),
            create_graph=True, retain_graph=True
        )[0]
        
        gradient_penalty = ((gradients.norm(2, dim=1) - 1) ** 2).mean()
        wasserstein_loss = torch.mean(self.discriminator(fake_data)) - torch.mean(self.discriminator(real_data))
        
        return wasserstein_loss + 10 * gradient_penalty

3.3 边界生成与验证

class BoundaryGenerator:
    def __init__(self, wgan_model, resolution=100):
        self.model = wgan_model
        self.resolution = resolution
        
    def generate_boundary(self, feature_ranges):
        """生成安全域边界"""
        boundaries = []
        
        # 在特征空间网格采样
        for i in range(len(feature_ranges)):
            for j in range(i+1, len(feature_ranges)):
                boundary_2d = self._generate_2d_boundary(
                    feature_ranges, i, j)
                boundaries.append(boundary_2d)
                
        return boundaries
    
    def _generate_2d_boundary(self, feature_ranges, dim1, dim2):
        """在二维平面上生成边界"""
        x = np.linspace(feature_ranges[dim1][0], feature_ranges[dim1][1], self.resolution)
        y = np.linspace(feature_ranges[dim2][0], feature_ranges[dim2][1], self.resolution)
        
        boundary_points = []
        for x_val in x:
            for y_val in y:
                # 构造特征向量
                sample = np.zeros(len(feature_ranges))
                sample[dim1] = x_val
                sample[dim2] = y_val
                
                # 使用生成器判断边界点
                with torch.no_grad():
                    sample_tensor = torch.FloatTensor(sample).to(self.model.device)
                    security_score = self.model.discriminator(sample_tensor.unsqueeze(0))
                    
                    # 边界点判别条件
                    if abs(security_score.item()) < 0.1:  # 接近决策边界
                        boundary_points.append([x_val, y_val])
                        
        return np.array(boundary_points)

4. 与传统方法的对比分析

4.1 最小二乘法拟合边界

传统方法常用最小二乘法拟合安全边界，假设边界为简单几何形状：

from scipy.optimize import least_squares

class LeastSquaresBoundary:
    def __init__(self, degree=2):
        self.degree = degree
        
    def fit(self, boundary_points):
        """使用最小二乘法拟合边界曲线"""
        x_data = boundary_points[:, 0]
        y_data = boundary_points[:, 1]
        
        # 多项式拟合
        coeffs = np.polyfit(x_data, y_data, self.degree)
        self.poly = np.poly1d(coeffs)
        
    def predict(self, x_values):
        return self.poly(x_values)

4.2 性能对比指标

通过多个指标对比两种方法的性能：

def evaluate_methods(true_boundary, wgan_boundary, ls_boundary):
    metrics = {}
    
    # 1. 拟合误差
    wgan_error = np.mean(np.min(np.linalg.norm(
        true_boundary[:, np.newaxis] - wgan_boundary, axis=2), axis=1))
    ls_error = np.mean(np.min(np.linalg.norm(
        true_boundary[:, np.newaxis] - ls_boundary, axis=2), axis=1))
    
    metrics['fitting_error'] = {'WGAN-GP': wgan_error, 'LeastSquares': ls_error}
    
    # 2. 计算时间对比
    metrics['computation_time'] = {
        'WGAN-GP': '秒级', 
        'LeastSquares': '分钟级', 
        'Pointwise': '小时级'
    }
    
    # 3. 边界平滑度
    wgan_smoothness = calculate_smoothness(wgan_boundary)
    ls_smoothness = calculate_smoothness(ls_boundary)
    metrics['smoothness'] = {'WGAN-GP': wgan_smoothness, 'LeastSquares': ls_smoothness}
    
    return metrics

def calculate_smoothness(boundary):
    """计算边界平滑度"""
    if len(boundary) < 2:
        return 0
    derivatives = np.diff(boundary, axis=0)
    return np.mean(np.linalg.norm(derivatives, axis=1))

4.3 对比结果分析

通过实验对比，WGAN-GP方法在多个维度表现优异：

5. 关键参数调优与模型训练

5.1 WGAN-GP超参数设置

模型性能对超参数敏感，需要仔细调优：

class HyperparameterTuner:
    def __init__(self, data_dim):
        self.data_dim = data_dim
        
    def grid_search(self, train_data, param_grid):
        best_score = float('inf')
        best_params = {}
        
        for latent_dim in param_grid['latent_dim']:
            for lr in param_grid['learning_rate']:
                for gp_lambda in param_grid['gp_lambda']:
                    
                    model = WGAN_GP(
                        latent_dim=latent_dim, 
                        data_dim=self.data_dim
                    )
                    
                    # 训练模型
                    train_loss = self.train_model(model, train_data)
                    
                    # 评估模型
                    score = self.evaluate_model(model, train_data)
                    
                    if score < best_score:
                        best_score = score
                        best_params = {
                            'latent_dim': latent_dim,
                            'learning_rate': lr,
                            'gp_lambda': gp_lambda
                        }
        
        return best_params, best_score

# 推荐参数范围
param_grid = {
    'latent_dim': [50, 100, 200],
    'learning_rate': [0.0001, 0.0005, 0.001],
    'gp_lambda': [1, 10, 100]
}

5.2 训练过程监控

训练过程中需要监控关键指标，确保模型收敛：

def monitor_training(g_losses, d_losses, boundary_accuracies):
    """监控训练过程"""
    import matplotlib.pyplot as plt
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
    
    # 损失曲线
    ax1.plot(g_losses, label='Generator Loss')
    ax1.plot(d_losses, label='Discriminator Loss')
    ax1.set_xlabel('Epoch')
    ax1.set_ylabel('Loss')
    ax1.legend()
    
    # 边界精度
    ax2.plot(boundary_accuracies)
    ax2.set_xlabel('Epoch')
    ax2.set_ylabel('Boundary Accuracy')
    
    plt.tight_layout()
    plt.show()

6. 实际应用中的问题与解决方案

6.1 数据不足问题

电力系统安全数据往往有限，需要数据增强：

class DataAugmentation:
    def __init__(self, original_data):
        self.data = original_data
        
    def smote_augmentation(self, k_neighbors=5):
        """使用SMOTE算法进行过采样"""
        from sklearn.neighbors import NearestNeighbors
        
        augmented_data = []
        n_samples, n_features = self.data.shape
        
        for i in range(n_samples):
            # 找到k个最近邻
            nn = NearestNeighbors(n_neighbors=k_neighbors+1)
            nn.fit(self.data)
            distances, indices = nn.kneighbors(self.data[i:i+1])
            
            # 生成新样本
            for j in range(1, k_neighbors+1):
                neighbor_idx = indices[0][j]
                alpha = np.random.random()
                new_sample = self.data[i] + alpha * (self.data[neighbor_idx] - self.data[i])
                augmented_data.append(new_sample)
                
        return np.vstack([self.data, np.array(augmented_data)])

6.2 边界模糊问题

安全边界附近样本判别模糊，需要特殊处理：

def boundary_refinement(initial_boundary, discriminator, refinement_steps=10):
    """边界精细化处理"""
    refined_boundary = initial_boundary.copy()
    
    for step in range(refinement_steps):
        for i in range(len(refined_boundary)):
            point = refined_boundary[i]
            
            # 在点周围采样，寻找更精确的边界点
            perturbations = np.random.normal(0, 0.01, (100, point.shape[0]))
            candidate_points = point + perturbations
            
            # 使用判别器评估边界性
            with torch.no_grad():
                scores = discriminator(torch.FloatTensor(candidate_points))
                boundary_scores = torch.abs(scores)
                best_idx = torch.argmin(boundary_scores)
                
                refined_boundary[i] = candidate_points[best_idx]
                
    return refined_boundary

6.3 模型稳定性保障

生产环境中需要确保模型稳定性：

class ModelStability:
    def __init__(self, model, backup_model=None):
        self.model = model
        self.backup = backup_model
        self.stability_history = []
        
    def check_stability(self, test_data, threshold=0.1):
        """检查模型输出稳定性"""
        with torch.no_grad():
            predictions = []
            for _ in range(10):  # 多次预测检验稳定性
                pred = self.model.discriminator(test_data)
                predictions.append(pred.cpu().numpy())
            
            std_dev = np.std(predictions)
            self.stability_history.append(std_dev)
            
            return std_dev < threshold
    
    def fallback_to_backup(self, current_data):
        """主模型不稳定时切换到备用模型"""
        if not self.check_stability(current_data):
            print("主模型不稳定，切换到备用模型")
            return self.backup
        return self.model

7. 生产环境部署建议

7.1 性能优化策略

实际部署时需要优化推理速度：

class InferenceOptimizer:
    def __init__(self, model):
        self.model = model
        
    def quantize_model(self):
        """模型量化加速"""
        model_quantized = torch.quantization.quantize_dynamic(
            self.model, {nn.Linear}, dtype=torch.qint8
        )
        return model_quantized
    
    def optimize_batch_size(self, test_data, max_batch_size=1024):
        """寻找最优批处理大小"""
        best_time = float('inf')
        best_batch_size = 1
        
        for batch_size in [1, 8, 32, 64, 128, 256, 512, 1024]:
            if batch_size > len(test_data):
                break
                
            start_time = time.time()
            with torch.no_grad():
                for i in range(0, len(test_data), batch_size):
                    batch = test_data[i:i+batch_size]
                    _ = self.model.discriminator(batch)
                    
            elapsed = time.time() - start_time
            if elapsed < best_time:
                best_time = elapsed
                best_batch_size = batch_size
                
        return best_batch_size

7.2 监控与告警

生产环境需要完善的监控体系：

class ProductionMonitor:
    def __init__(self, model, acceptable_drift=0.05):
        self.model = model
        self.baseline_performance = None
        self.drift_threshold = acceptable_drift
        
    def establish_baseline(self, validation_data):
        """建立性能基线"""
        with torch.no_grad():
            predictions = self.model.discriminator(validation_data)
            self.baseline_performance = {
                'mean_score': torch.mean(predictions).item(),
                'std_score': torch.std(predictions).item()
            }
    
    def check_concept_drift(self, new_data):
        """检查概念漂移"""
        with torch.no_grad():
            new_predictions = self.model.discriminator(new_data)
            new_mean = torch.mean(new_predictions).item()
            
            drift_amount = abs(new_mean - self.baseline_performance['mean_score'])
            return drift_amount > self.drift_threshold

WGAN-GP方法为动态安全域边界生成提供了新的技术路径，相比传统方法在计算效率和边界质量方面都有显著提升。实际应用中需要根据具体电力系统特点调整特征选择和模型参数，同时建立完善的监控机制确保在线应用的可靠性。