# 分布式鲁棒优化实战：用Python处理风电预测误差（附KL散度代码） 在能源行业，尤其是风电、光伏等可再生能源大规模并网的今天，预测误差带来的不确定性已经成为调度、规划和交易中无法回避的核心挑战。传统的随机规划方法依赖于对不确定性参数（如风电出力）概率分布的精确假设，这在现实中往往难以满足——历史数据有限，未来分布未知，一个“错误”的分布假设可能导致优化结果在实际运行中失效。而经典的鲁棒优化虽然对分布信息要求极低，但其“最坏情况”的保守性又常常以牺牲巨大的经济性为代价。 有没有一种方法，能在“已知部分信息”和“避免过度保守”之间找到平衡点？这正是**分布式鲁棒优化**（Distributionally Robust Optimization, DRO）大显身手的领域。它不假设一个精确的分布，而是构建一个包含所有“可能”分布的**模糊集**（Ambiguity Set），并在这个集合中优化最坏情况下的期望性能。对于风电这类历史数据相对丰富、但未来分布难以精确刻画的不确定性，DRO提供了一种既利用数据、又承认数据局限性的务实框架。 本文面向能源行业的从业者、分析师和研究者，旨在抛开复杂的理论推导，通过Python代码实战，带你走完一个完整的DRO应用流程：从**历史风电数据清洗**，到基于**KL散度**构建模糊集，再到构建并求解一个考虑风电预测误差的**分布鲁棒优化模型**。我们将使用`pandas`进行数据处理，用`scipy`进行优化求解，并附上可直接运行的代码块。你会发现，DRO并非遥不可及的理论，而是可以快速上手的强大工具。 ## 1. 问题场景与数据准备：风电预测误差的挑战 假设我们负责一个包含风电场的电力系统日前经济调度。每天需要基于风电功率预测来决定常规机组的启停和出力计划。风电预测存在误差，如果计划过于乐观（低估了实际风电），可能导致功率缺额，需要调用昂贵的备用或造成切负荷；如果计划过于保守（高估了实际风电），则可能让煤电机组多发电，造成不必要的燃料成本和碳排放。 我们拥有过去一段时间（例如30天）的风电实际出力与预测出力的历史数据。**核心问题是：如何利用这有限的历史数据，构建一个关于预测误差的“不确定性模型”，使得做出的调度决策既能抵御合理的误差波动，又不会为极端小概率事件付出过高成本？** 首先，我们模拟生成一段历史数据，并观察预测误差的特性。 ```python import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from scipy import stats, optimize import warnings warnings.filterwarnings('ignore') # 设置随机种子以保证结果可复现 np.random.seed(42) # 模拟30天，每15分钟一个点（每天96个点）的风电预测和实际值 n_days = 30 points_per_day = 96 n_points = n_days * points_per_day # 生成一个简单的日周期波动作为预测值（假设预测是完美的周期模式） time_index = pd.date_range(start='2023-01-01', periods=n_points, freq='15T') daily_pattern = 50 + 30 * np.sin(2 * np.pi * np.arange(points_per_day) / points_per_day) # 日均值50MW，振幅30MW wind_forecast = np.tile(daily_pattern, n_days) # 生成实际值：在预测值基础上加上随机误差。误差分布我们假设是均值为0，但真实分布未知。 # 为了模拟DRO的价值，我们让误差分布不是完美的正态分布，而是带有一定的偏态和厚尾。 shape_param = 5 # 形状参数，控制偏态 scale_param = 15 # 尺度参数，控制误差大小 # 使用Beta分布平移缩放来生成非对称误差 base_error = np.random.beta(shape_param, shape_param, n_points) - 0.5 # 均值为0的Beta分布 wind_error = base_error * scale_param * (1 + 0.3 * np.random.randn(n_points)) # 加入随机波动 wind_actual = wind_forecast + wind_error # 确保实际出力非负 wind_actual = np.maximum(wind_actual, 0) # 创建DataFrame df_wind = pd.DataFrame({ 'time': time_index, 'forecast_MW': wind_forecast, 'actual_MW': wind_actual, 'error_MW': wind_actual - wind_forecast }) df_wind['day'] = df_wind['time'].dt.date print("风电数据前5行：") print(df_wind.head()) print(f"\n数据统计描述（误差）：") print(df_wind['error_MW'].describe()) # 可视化第一天的预测与实际值 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) day_one = df_wind[df_wind['day'] == df_wind['day'].unique()[0]] axes[0, 0].plot(day_one['time'], day_one['forecast_MW'], label='预测', linewidth=2) axes[0, 0].plot(day_one['time'], day_one['actual_MW'], label='实际', alpha=0.7) axes[0, 0].set_title('第一天风电出力预测 vs 实际') axes[0, 0].set_ylabel('功率 (MW)') axes[0, 0].legend() axes[0, 0].grid(True, linestyle='--', alpha=0.7) # 误差分布直方图 axes[0, 1].hist(df_wind['error_MW'], bins=50, density=True, alpha=0.6, color='skyblue', edgecolor='black') axes[0, 1].axvline(x=0, color='red', linestyle='--', label='误差均值') axes[0, 1].set_title('风电预测误差分布直方图') axes[0, 1].set_xlabel('误差 (MW)') axes[0, 1].set_ylabel('密度') axes[0, 1].legend() axes[0, 1].grid(True, linestyle='--', alpha=0.7) # 误差的Q-Q图（检验正态性） stats.probplot(df_wind['error_MW'], dist="norm", plot=axes[1, 0]) axes[1, 0].set_title('误差Q-Q图 (vs 正态分布)') axes[1, 0].grid(True, linestyle='--', alpha=0.7) # 误差自相关图（检查时间相关性） from pandas.plotting import autocorrelation_plot autocorrelation_plot(df_wind['error_MW'], ax=axes[1, 1]) axes[1, 1].set_title('误差自相关图') axes[1, 1].set_ylim([-0.2, 0.2]) axes[1, 1].grid(True, linestyle='--', alpha=0.7) plt.tight_layout() plt.show() ``` 运行以上代码，我们会得到风电预测误差的初步分析。从直方图和Q-Q图可以明显看出，误差分布并非完美的正态分布，存在一定的偏态。自相关图则显示误差在短时间内可能存在微弱的相关性，但为了简化模型，我们通常假设各时刻误差独立。**正是这种分布特性的不确定性，使得我们无法放心地使用基于正态假设的随机规划，而DRO的模糊集可以自然地容纳这种分布的不确定性。** > **提示**：在实际项目中，你需要用真实的历史数据替换这里的模拟数据。数据清洗步骤可能包括处理缺失值、剔除异常值（如由于风机故障导致的长时间零出力）、以及可能需要的误差归一化（如将误差表示为预测值的百分比）。 ## 2. 模糊集构建：从经验分布到KL散度球 分布式鲁棒优化的核心在于**模糊集**的构建。模糊集定义了我们认为“可能”的真实概率分布集合。构建方式主要分两类：基于矩信息（如均值、方差已知，但具体分布未知）和基于统计距离（如与经验分布的距离不超过某个阈值）。后者在大数据背景下更为常用，因为它更直接地利用了历史数据。 **KL散度**（Kullback-Leibler Divergence），又称相对熵，是衡量两个概率分布差异的常用工具。给定经验分布 $\hat{\mathbb{P}}$（由N个历史误差样本构成，每个样本权重为1/N）和一个待考察的分布 $\mathbb{P}$，它们之间的KL散度定义为： $$ D_{KL}(\mathbb{P} || \hat{\mathbb{P}}) = \sum_{\omega=1}^{N} p_\omega \log\left(\frac{p_\omega}{1/N}\right) $$ 其中 $p_\omega$ 是分布 $\mathbb{P}$ 赋予第 $\omega$ 个历史样本的概率。KL散度非负，且当且仅当 $\mathbb{P} = \hat{\mathbb{P}}$ 时为零。 基于KL散度的模糊集构造如下： $$ \mathcal{F}_{KL}(\theta) = \left\{ \mathbb{P} \in \mathcal{P}(\Xi) : D_{KL}(\mathbb{P} || \hat{\mathbb{P}}) \leq \theta \right\} $$ 其中 $\theta \geq 0$ 是**半径**或**容忍度**参数。$\theta$ 越大，模糊集包含的分布越多，模型越保守；$\theta = 0$ 时，模糊集退化为只包含经验分布，DRO退化为随机规划（基于样本平均近似）。 下面，我们用代码计算经验分布，并可视化KL散度模糊集的含义。 ```python # 提取历史误差样本 historical_errors = df_wind['error_MW'].values N = len(historical_errors) print(f"历史样本数量 N = {N}") # 经验分布：每个样本权重为 1/N empirical_probs = np.ones(N) / N # 定义一个函数，计算离散分布q相对于经验分布p_emp的KL散度 def kl_divergence(q, p_emp): """ 计算离散概率向量q相对于p_emp的KL散度。 假设q和p_emp都是归一化的概率向量，且维度相同。 处理q中可能为0的元素（定义0*log(0)=0）。 """ # 避免log(0)的情况 mask = q > 0 return np.sum(q[mask] * np.log(q[mask] / p_emp[mask])) # 示例：生成一个扰动后的分布，计算其KL散度 np.random.seed(123) # 对经验概率进行随机扰动，然后重新归一化 perturbation = 0.1 * np.random.randn(N) q_test = empirical_probs * np.exp(perturbation) q_test = q_test / q_test.sum() # 归一化 kl_value = kl_divergence(q_test, empirical_probs) print(f"一个扰动分布的KL散度值: {kl_value:.6f}") # 探索不同扰动强度下的KL散度 def generate_perturbed_distribution(base_probs, perturbation_scale): """生成一个扰动分布，扰动尺度由perturbation_scale控制。""" perturbation = perturbation_scale * np.random.randn(N) # 使用softmax形式的扰动以保证正性和归一化（更稳定） log_probs = np.log(base_probs + 1e-10) + perturbation q = np.exp(log_probs) return q / q.sum() # 采样多个扰动分布，计算其KL散度 np.random.seed(456) theta_values = [] perturbation_scales = [0.01, 0.05, 0.1, 0.2, 0.5] for scale in perturbation_scales: for _ in range(20): # 每个尺度生成20个样本 q_pert = generate_perturbed_distribution(empirical_probs, scale) theta_values.append((scale, kl_divergence(q_pert, empirical_probs))) theta_df = pd.DataFrame(theta_values, columns=['perturbation_scale', 'KL_divergence']) # 可视化 plt.figure(figsize=(10, 6)) for scale in perturbation_scales: subset = theta_df[theta_df['perturbation_scale'] == scale]['KL_divergence'] plt.scatter([scale]*len(subset), subset, alpha=0.6, label=f'scale={scale}') plt.axhline(y=0, color='k', linestyle='--', alpha=0.3) plt.xlabel('扰动尺度 (Perturbation Scale)') plt.ylabel('KL散度 $D_{KL}$') plt.title('不同扰动强度下生成分布相对于经验分布的KL散度') plt.legend() plt.grid(True, linestyle='--', alpha=0.7) plt.show() # 选择一个合理的theta值。一种常见方法是基于样本量进行设置。 # 例如，根据大数定律，经验分布接近真实分布的概率与样本量和theta有关。 # 这里我们采用一个经验性选择：theta = 0.1 / N 或类似。我们也可以使用交叉验证来确定。 theta_selected = 0.05 # 这是一个需要根据风险偏好调整的超参数 print(f"\n选择的KL散度模糊集半径 theta = {theta_selected}") ``` 通过上述可视化，我们可以直观感受扰动尺度与KL散度大小的关系。**半径 $\theta$ 的选择是DRO应用中的关键**，它体现了决策者的风险厌恶程度。$\theta$ 越大，考虑的分布不确定性越广，模型越保守，成本通常越高。在实际应用中，可以通过**交叉验证**或**回溯测试**，在历史数据上测试不同 $\theta$ 值对应的决策表现，从而选择一个在鲁棒性和经济性之间平衡的值。 ## 3. 构建分布鲁棒优化模型：以简单经济调度为例 现在，我们将风电预测误差的不确定性，通过KL散度模糊集，嵌入到一个简化的电力系统经济调度模型中。考虑一个单时段问题：系统总负荷为 $L$，有一台常规火电机组，其发电成本为二次函数 $C(p_g) = a p_g^2 + b p_g + c$，出力范围 $[p_g^{min}, p_g^{max}]$。风电预测值为 $w_f$，但实际出力 $w$ 不确定，其分布属于模糊集 $\mathcal{F}_{KL}(\theta)$。我们需要决定火电机组出力 $p_g$，以满足负荷平衡，并最小化**最坏情况下的期望总成本**。 负荷平衡约束为：$p_g + w = L$。由于 $w$ 不确定，这是一个随机约束。在DRO框架下，我们要求该约束对于模糊集 $\mathcal{F}_{KL}(\theta)$ 中的所有分布**几乎必然成立**（即鲁棒约束），或者以一定概率成立（机会约束）。这里我们先考虑**鲁棒约束**，即对于模糊集中所有可能的分布，约束都必须满足。这意味着我们需要针对最坏情况下的风电出力（即可能的最小值）来安排火电出力，但“最坏情况”是在模糊集定义的分布范围内寻找。 为了简化并突出DRO的核心，我们考虑一个更易处理的**期望成本最小化**问题，并将不确定性体现在目标函数中： $$ \min_{p_g} \sup_{\mathbb{P} \in \mathcal{F}_{KL}(\theta)} \mathbb{E}_{\mathbb{P}}[C(p_g) + \lambda \cdot \max(L - p_g - w, 0)] $$ 其中第二项是惩罚可能的功率缺额（失负荷），$\lambda$ 是失负荷惩罚系数（远大于发电成本系数）。这个问题等价于： $$ \min_{p_g, t} \sup_{\mathbb{P} \in \mathcal{F}_{KL}(\theta)} \mathbb{E}_{\mathbb{P}}[C(p_g) + \lambda \cdot s] $$ 约束为：$s \geq L - p_g - w$，且 $s \geq 0$。 **DRO问题的求解技巧在于利用对偶理论**。对于基于KL散度的模糊集，上述内层“sup”问题可以转化为一个**凸优化问题**，并且有解析解或易于求解的形式。具体来说，有下面的重要结论： > **定理**：对于基于KL散度的模糊集 $\mathcal{F}_{KL}(\theta)$，优化问题 $\sup_{\mathbb{P} \in \mathcal{F}_{KL}(\theta)} \mathbb{E}_{\mathbb{P}}[f(w)]$ 的最优值等于： > $$ > \inf_{\beta \geq 0} \left\{ \beta \theta + \beta \log \left( \frac{1}{N} \sum_{\omega=1}^{N} \exp(f(w_\omega) / \beta) \right) \right\} > $$ > 其中 $w_\omega$ 是历史样本。 利用这个定理，我们可以将原问题转化为一个关于 $p_g, s, \beta$ 的确定性优化问题。下面我们用代码实现这个转化和求解过程。 ```python # 定义系统参数 L = 200.0 # 系统总负荷，单位MW a, b, c = 0.05, 20, 100 # 火电机组成本系数 C(p)=a*p^2 + b*p + c p_min, p_max = 50, 250 # 火电机组出力上下限，单位MW lambda_voLL = 1000 # 失负荷价值 (Value of Lost Load)，单位 $/MWh，远高于发电成本 # 风电预测值（当前时段） w_forecast_current = 60.0 # MW # 注意：我们的历史误差样本是针对整个时间序列的。在单时段问题中，我们假设当前时段的误差分布与历史同期或整体类似。 # 这里我们直接使用全部历史误差样本来构建模糊集，相当于假设误差分布是平稳的。 # 更精细的做法是选取相似时段（如相同小时、相似天气）的历史误差。 # 为当前时段构建风电出力的可能场景：预测值 + 历史误差 w_scenarios = w_forecast_current + historical_errors # 检查是否有异常场景（如风电出力为负），将其修正为0 w_scenarios = np.maximum(w_scenarios, 0) N = len(w_scenarios) print(f"场景数量: {N}") print(f"风电场景范围: [{w_scenarios.min():.2f}, {w_scenarios.max():.2f}] MW") # 定义基于KL散度DRO的优化问题求解函数 def solve_dro_kl(theta, w_scenarios, lambda_voLL, a, b, c, L, p_min, p_max): """ 求解单时段分布鲁棒经济调度问题。 目标：最小化最坏情况期望成本，包括火电成本和失负荷惩罚。 模糊集：基于KL散度，半径theta。 返回：最优火电出力 p_g_opt, 最优目标值 obj_opt, 最优对偶变量 beta_opt """ N = len(w_scenarios) # 定义内层函数：给定p_g和beta，计算目标值（利用定理） def objective_for_beta(x): p_g, beta = x[0], x[1] if p_g < p_min or p_g > p_max or beta <= 0: return np.inf # 计算每个场景下的总成本 f(w) = C(p_g) + lambda * max(L - p_g - w, 0) # 注意：C(p_g)是确定项，可以提到期望外面。但为了公式统一，我们放在里面。 cost_deterministic = a * p_g**2 + b * p_g + c # 计算每个场景的失负荷量 load_shed = np.maximum(L - p_g - w_scenarios, 0) f_vals = cost_deterministic + lambda_voLL * load_shed # 当beta很小时，指数项可能爆炸，需要数值处理 # 计算 log-sum-exp 项 max_f = np.max(f_vals) exp_terms = np.exp((f_vals - max_f) / beta) # 减去最大值以提高数值稳定性 log_sum_exp = max_f + beta * np.log(np.mean(exp_terms)) # 最终目标值：beta * theta + beta * log( (1/N) sum exp(f/beta) ) obj_val = beta * theta + log_sum_exp return obj_val # 使用优化器求解关于p_g和beta的最小化问题 # 初始猜测 p_g_init = np.clip(L - w_forecast_current, p_min, p_max) # 基于预测的平衡点 beta_init = 10.0 # 初始对偶变量 bounds = [(p_min, p_max), (1e-6, 1000)] # beta必须为正 result = optimize.minimize(objective_for_beta, x0=[p_g_init, beta_init], bounds=bounds, method='L-BFGS-B', options={'ftol': 1e-9, 'maxiter': 1000}) if not result.success: print(f"优化警告: {result.message}") p_g_opt, beta_opt = result.x obj_opt = result.fun return p_g_opt, beta_opt, obj_opt # 定义作为对比的基准方法求解函数 def solve_deterministic(w_realized): """确定性优化：假设知道风电实现值""" p_g = np.clip(L - w_realized, p_min, p_max) cost = a * p_g**2 + b * p_g + c + lambda_voLL * max(L - p_g - w_realized, 0) return p_g, cost def solve_stochastic_erm(w_scenarios): """随机规划（期望值模型）：最小化平均成本""" # 这是一个两阶段随机规划，但这里简化为直接对场景平均求最优 # 通过枚举p_g，计算平均成本 p_grid = np.linspace(p_min, p_max, 200) avg_costs = [] for p_g in p_grid: load_shed = np.maximum(L - p_g - w_scenarios, 0) total_costs = a * p_g**2 + b * p_g + c + lambda_voLL * load_shed avg_costs.append(np.mean(total_costs)) idx_opt = np.argmin(avg_costs) p_g_opt = p_grid[idx_opt] obj_opt = avg_costs[idx_opt] return p_g_opt, obj_opt def solve_robust_wc(w_scenarios): """经典鲁棒优化：针对最坏场景（最小风电出力）""" w_worst = w_scenarios.min() p_g_opt = np.clip(L - w_worst, p_min, p_max) # 计算在最坏场景下的成本（虽然决策是针对最坏场景的，但成本计算时我们仍用该场景） cost = a * p_g_opt**2 + b * p_g_opt + c + lambda_voLL * max(L - p_g_opt - w_worst, 0) return p_g_opt, cost # 求解不同方法 theta = 0.05 p_g_dro, beta_dro, obj_dro = solve_dro_kl(theta, w_scenarios, lambda_voLL, a, b, c, L, p_min, p_max) p_g_erm, obj_erm = solve_stochastic_erm(w_scenarios) p_g_wc, obj_wc = solve_robust_wc(w_scenarios) print("=== 不同优化方法结果比较 ===") print(f"风电预测值: {w_forecast_current:.2f} MW") print(f"系统负荷: {L:.2f} MW") print("-" * 50) print(f"分布鲁棒优化 (DRO, theta={theta}):") print(f" 火电最优出力: {p_g_dro:.2f} MW") print(f" 最坏情况期望成本: {obj_dro:.2f} $") print(f" 对偶变量 beta: {beta_dro:.6f}") print() print(f"随机规划 (期望值模型, ERM):") print(f" 火电最优出力: {p_g_erm:.2f} MW") print(f" 平均成本: {obj_erm:.2f} $") print() print(f"经典鲁棒优化 (最坏场景, WC):") print(f" 火电最优出力: {p_g_wc:.2f} MW") print(f" 最坏场景成本: {obj_wc:.2f} $") print() # 计算一些关键指标 # 计算DRO决策在不同场景下的实际成本 load_shed_dro = np.maximum(L - p_g_dro - w_scenarios, 0) costs_dro = a * p_g_dro**2 + b * p_g_dro + c + lambda_voLL * load_shed_dro avg_cost_dro = np.mean(costs_dro) worst_cost_dro = np.max(costs_dro) # 计算ERM决策在不同场景下的实际成本 load_shed_erm = np.maximum(L - p_g_erm - w_scenarios, 0) costs_erm = a * p_g_erm**2 + b * p_g_erm + c + lambda_voLL * load_shed_erm avg_cost_erm = np.mean(costs_erm) worst_cost_erm = np.max(costs_erm) print("=== 决策在不同场景下的表现（样本外）===") print(f" DRO决策 ERM决策") print(f"平均成本 {avg_cost_dro:.2f} $ {avg_cost_erm:.2f} $") print(f"最坏成本 {worst_cost_dro:.2f} $ {worst_cost_erm:.2f} $") print(f"失负荷概率 {(load_shed_dro > 0).mean():.2%} {(load_shed_erm > 0).mean():.2%}") print(f"平均失负荷量 {load_shed_dro.mean():.2f} MW {load_shed_erm.mean():.2f} MW") ``` 运行上述代码，我们可以立即比较DRO、随机规划（ERM）和经典鲁棒优化（WC）三种方法的结果。通常你会发现： * **ERM**追求平均成本最低，但可能在少数恶劣场景下（风电出力极低）产生高额的失负荷成本，风险较高。 * **WC**完全规避风险，但火电安排过高，导致在大多数场景下发电成本偏高，经济性差。 * **DRO**通过调节 $\theta$，在两者之间取得平衡。当 $\theta$ 较小时，结果接近ERM；$\theta$ 较大时，结果接近WC。 ## 4. 深入分析：theta的影响与模型扩展 为了更全面地理解DRO，我们需要系统地研究模糊集半径 $\theta$ 对决策和性能的影响。同时，我们将模型扩展到一个更实用的**多时段**场景，并引入**机会约束**。 ### 4.1 分析theta的敏感性 $\theta$ 是DRO模型的“调节旋钮”。我们通过改变 $\theta$，观察最优决策、成本以及风险指标的变化。 ```python # 研究theta对DRO结果的影响 theta_range = np.logspace(-3, 1, 20) # 从0.001到10，对数间隔 results = [] for theta_val in theta_range: p_g_opt, beta_opt, obj_opt = solve_dro_kl(theta_val, w_scenarios, lambda_voLL, a, b, c, L, p_min, p_max) # 计算该决策下的样本平均成本和最坏成本 load_shed = np.maximum(L - p_g_opt - w_scenarios, 0) costs = a * p_g_opt**2 + b * p_g_opt + c + lambda_voLL * load_shed avg_cost = np.mean(costs) worst_cost = np.max(costs) prob_loss = (load_shed > 0).mean() results.append({ 'theta': theta_val, 'p_g_opt': p_g_opt, 'beta_opt': beta_opt, 'obj_dro': obj_opt, 'avg_cost': avg_cost, 'worst_cost': worst_cost, 'prob_loss': prob_loss }) results_df = pd.DataFrame(results) # 可视化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) ax = axes[0, 0] ax.semilogx(results_df['theta'], results_df['p_g_opt'], 'o-', linewidth=2, markersize=6) ax.axhline(y=p_g_erm, color='r', linestyle='--', label=f'ERM决策 ({p_g_erm:.1f} MW)') ax.axhline(y=p_g_wc, color='g', linestyle='--', label=f'WC决策 ({p_g_wc:.1f} MW)') ax.set_xlabel(r'模糊集半径 $\theta$') ax.set_ylabel('火电最优出力 (MW)') ax.set_title(r'$\theta$ 对最优决策的影响') ax.legend() ax.grid(True, which="both", ls="--", alpha=0.5) ax = axes[0, 1] ax.semilogx(results_df['theta'], results_df['avg_cost'], 's-', label='平均成本', linewidth=2) ax.semilogx(results_df['theta'], results_df['worst_cost'], '^-', label='最坏成本', linewidth=2) ax.axhline(y=obj_erm, color='r', linestyle='--', label=f'ERM平均成本 ({obj_erm:.1f} $)') ax.axhline(y=obj_wc, color='g', linestyle='--', label=f'WC最坏成本 ({obj_wc:.1f} $)') ax.set_xlabel(r'模糊集半径 $\theta$') ax.set_ylabel('成本 ($)') ax.set_title(r'$\theta$ 对成本的影响') ax.legend() ax.grid(True, which="both", ls="--", alpha=0.5) ax = axes[1, 0] ax.semilogx(results_df['theta'], results_df['beta_opt'], 'd-', linewidth=2) ax.set_xlabel(r'模糊集半径 $\theta$') ax.set_ylabel(r'对偶变量 $\beta$') ax.set_title(r'$\theta$ 对最优对偶变量 $\beta$ 的影响') ax.grid(True, which="both", ls="--", alpha=0.5) ax = axes[1, 1] ax.semilogx(results_df['theta'], results_df['prob_loss'], 'v-', linewidth=2, color='purple') ax.set_xlabel(r'模糊集半径 $\theta$') ax.set_ylabel('失负荷概率') ax.set_title(r'$\theta$ 对失负荷风险的影响') ax.grid(True, which="both", ls="--", alpha=0.5) plt.tight_layout() plt.show() # 输出关键观察 print("关键观察：") print("1. 当 theta -> 0 时，DRO决策趋近于ERM决策，平均成本最低，但最坏成本和失负荷风险较高。") print("2. 当 theta 增大时，火电出力增加，平均成本上升，但最坏成本和失负荷风险下降。") print("3. 当 theta -> 无穷大时，DRO决策趋近于WC决策，完全规避风险，但经济性最差。") print("4. 对偶变量 beta 随着 theta 增大而增大，反映了对不确定性的‘定价’。") ``` ### 4.2 扩展到多时段与机会约束 单时段鲁棒约束可能过于严格。在实际调度中，我们允许一定的失负荷风险，这便引入了**机会约束**。分布鲁棒机会约束（DRCC）要求：在最坏的可能分布下，约束被满足的概率不低于一个置信水平 $1-\epsilon$。 对于基于KL散度的模糊集，一个幸运的结论是：分布鲁棒机会约束可以等价转化为一个**修正了置信水平**的传统机会约束（基于经验分布）。具体来说，约束 $$ \inf_{\mathbb{P} \in \mathcal{F}_{KL}(\theta)} \mathbb{P}(p_g + w \geq L) \geq 1 - \epsilon $$ 等价于 $$ \hat{\mathbb{P}}(p_g + w \geq L) \geq 1 - \bar{\epsilon} $$ 其中 $\bar{\epsilon} = 1 - \inf_{t \in (0,1)} \frac{e^{-\theta} t^{1-\epsilon} - 1}{t - 1}$。这意味着我们只需要对经验分布施加一个更严格的概率约束即可。 下面，我们实现一个考虑机会约束的多时段简化模型。假设有T个时段，每个时段的风电预测误差独立，且拥有各自的历史误差样本（这里为简化，使用同一组样本）。目标是最小化总期望成本，同时满足每个时段负荷平衡的分布鲁棒机会约束。 ```python def solve_multiperiod_drcc(theta, epsilon, w_forecast_series, historical_errors_multi, lambda_voLL, a, b, c, L_series, p_min, p_max): """ 求解多时段分布鲁棒机会约束经济调度。 假设各时段独立，使用相同的theta和epsilon。 w_forecast_series: 各时段风电预测值列表，长度T。 historical_errors_multi: 列表，每个元素是一个数组，包含该时段的历史误差样本。 L_series: 各时段负荷列表，长度T。 返回: 各时段最优火电出力列表 p_g_opt_series, 总成本 """ T = len(w_forecast_series) # 计算修正后的置信水平 epsilon_bar # 根据公式 epsilon_bar = 1 - inf_{t in (0,1)} [ (e^{-theta} * t^{1-epsilon} - 1) / (t - 1) ] # 可以通过一维优化求解 def func(t): if t <= 0 or t >= 1: return np.inf numerator = np.exp(-theta) * (t ** (1 - epsilon)) - 1 denominator = t - 1 return 1 - numerator / denominator t_vals = np.linspace(1e-6, 1-1e-6, 1000) epsilon_bar_vals = np.array([func(t) for t in t_vals]) epsilon_bar = epsilon_bar_vals.min() print(f"原始置信水平 1-epsilon = {1-epsilon:.3f}, 修正后置信水平 1-epsilon_bar = {1-epsilon_bar:.3f}") p_g_opt_series = [] total_cost = 0 for t in range(T): w_forecast = w_forecast_series[t] errors = historical_errors_multi[t] # 该时段的历史误差样本 L = L_series[t] w_scenarios = w_forecast + errors w_scenarios = np.maximum(w_scenarios, 0) N_t = len(w_scenarios) # 机会约束等价于：经验分布下，p_g + w >= L 的概率 >= 1 - epsilon_bar # 即：在历史样本中，满足 p_g + w_i >= L 的比例 >= 1 - epsilon_bar # 这等价于：p_g >= L - w_{(k)}，其中 w_{(k)} 是风电场景的 (1-epsilon_bar) 分位数（从小到大排序） # 注意：这里是针对风电出力的下分位数，因为我们需要 p_g + w >= L，即 p_g >= L - w。 # 所以，我们需要风电的 (epsilon_bar) 分位数（因为小风不利）。 w_sorted = np.sort(w_scenarios) index_k = int(np.floor(epsilon_bar * N_t)) # 下分位数索引 if index_k >= N_t: index_k = N_t - 1 w_epsilon_bar = w_sorted[index_k] # 满足机会约束的最小火电出力 p_g_min_required = L - w_epsilon_bar # 同时还要满足技术出力上下限 p_g_opt = np.clip(p_g_min_required, p_min, p_max) # 计算该时段的期望成本（基于经验分布） load_shed = np.maximum(L - p_g_opt - w_scenarios, 0) costs = a * p_g_opt**2 + b * p_g_opt + c + lambda_voLL * load_shed period_cost = np.mean(costs) p_g_opt_series.append(p_g_opt) total_cost += period_cost # 打印该时段信息 print(f"时段 {t+1}: 预测风电={w_forecast:.1f} MW, 负荷={L:.1f} MW, " f"风电{epsilon_bar:.1%}分位数={w_epsilon_bar:.1f} MW, " f"所需最小火电={p_g_min_required:.1f} MW, 实际安排火电={p_g_opt:.1f} MW, " f"时段期望成本={period_cost:.1f} $") return p_g_opt_series, total_cost # 模拟多时段数据 T = 24 # 24小时 np.random.seed(888) # 生成负荷曲线（日模式） load_base = 180 load_amp = 40 L_series = load_base + load_amp * np.sin(2 * np.pi * np.arange(T) / T + np.pi/4) # 生成风电预测曲线 wind_base = 50 wind_amp = 30 w_forecast_series = wind_base + wind_amp * np.sin(2 * np.pi * np.arange(T) / T - np.pi/6) # 为每个时段生成独立的历史误差样本（这里简单地从总误差中随机抽取） historical_errors_multi = [] for t in range(T): # 假设每个时段有100个历史误差样本 indices = np.random.choice(len(historical_errors), size=100, replace=True) historical_errors_multi.append(historical_errors[indices]) # 设置机会约束参数 epsilon = 0.05 # 允许5%的失负荷风险 theta_drcc = 0.1 print("=== 多时段分布鲁棒机会约束调度 ===") print(f"参数: theta={theta_drcc}, 允许失负荷风险 epsilon={epsilon}") print("-" * 80) p_g_opt_drcc, total_cost_drcc = solve_multiperiod_drcc(theta_drcc, epsilon, w_forecast_series, historical_errors_multi, lambda_voLL, a, b, c, L_series, p_min, p_max) print("-" * 80) print(f"多时段总期望成本: {total_cost_drcc:.2f} $") # 可视化多时段调度计划 fig, axes = plt.subplots(2, 1, figsize=(14, 10)) t_range = np.arange(1, T+1) axes[0].plot(t_range, L_series, 'k-', linewidth=3, label='系统负荷') axes[0].plot(t_range, w_forecast_series, 'b-', linewidth=2, label='风电预测', alpha=0.7) axes[0].plot(t_range, p_g_opt_drcc, 'r-', linewidth=2, label='火电计划 (DRCC)') axes[0].fill_between(t_range, 0, p_g_opt_drcc, color='red', alpha=0.2) axes[0].set_xlabel('时段 (小时)') axes[0].set_ylabel('功率 (MW)') axes[0].set_title('多时段分布鲁棒机会约束调度计划') axes[0].legend() axes[0].grid(True, linestyle='--', alpha=0.7) # 计算净负荷（负荷-风电预测）和火电计划的差值 net_load = L_series - w_forecast_series axes[1].bar(t_range, net_load, width=0.6, label='净负荷 (负荷-风电预测)', color='gray', alpha=0.5) axes[1].plot(t_range, p_g_opt_drcc, 'ro-', linewidth=2, markersize=6, label='火电计划 (DRCC)') axes[1].set_xlabel('时段 (小时)') axes[1].set_ylabel('功率 (MW)') axes[1].set_title('净负荷与火电计划对比') axes[1].legend() axes[1].grid(True, linestyle='--', alpha=0.7) plt.tight_layout() plt.show() ``` 通过多时段DRCC模型，我们得到了一个考虑风险约束的调度计划。火电出力在净负荷较高的时段（如傍晚负荷高峰且风电下降时）被适当抬高，以抵御风电预测误差可能带来的风险，同时在风电充足的时段降低出力以节省成本。**修正后的置信水平 $1-\bar{\epsilon}$ 比原始要求 $1-\epsilon$ 更严格**，这正是为应对分布不确定性所付出的“鲁棒性溢价”。 ## 5. 工程实践要点与高级话题 在实际项目中应用分布鲁棒优化，以下几个要点需要特别关注： 1. **数据质量与代表性**：历史误差样本的质量直接决定模糊集的有效性。需要确保样本是独立同分布的，或至少是平稳的。对于风电，通常需要考虑误差与预测值、季节、天气等因素的相关性，可能需要构建条件模糊集。 2. **theta的校准**：半径 $\theta$ 是核心超参数。除了交叉验证，还可以使用**Wasserstein距离**等构建具有有限样本保证的模糊集，其半径 $\theta$ 可以根据样本量 $N$ 和期望的置信水平 $\eta$ 通过公式 $\theta \propto N^{-1/d}$ 来确定（$d$ 为随机变量维度），这提供了更强的统计解释。 3. **计算效率**：对于大规模问题（如多节点电网、多类型不确定性），直接使用基于KL散度的对偶变换可能导致非线性问题。此时，基于**Wasserstein距离**的DRO虽然模糊集构造不同，但常可转化为线性规划或二阶锥规划，更适合大规模求解。此外，**场景削减技术**（如K-means聚类）可用于减少历史样本数量 $N$，从而降低计算负担。 4. **与其他不确定性的结合**：电力系统中还存在其他不确定性，如负荷预测误差、设备故障、市场价格波动等。DRO可以处理多重不确定性，但模糊集的构造会变得更加复杂（如构建矩模糊集处理相关的不确定性）。 为了更直观地对比不同模糊集构建方法，下表总结了基于矩、KL散度和Wasserstein距离的三种典型DRO方法的特点： | 特性 | 基于矩的DRO | 基于KL散度的DRO | 基于Wasserstein距离的DRO | | :--- | :--- | :--- | :--- | | **所需信息** | 随机变量的矩（如均值、方差）的不确定范围 | 历史数据样本（经验分布）及散度容差 $\theta$ | 历史数据样本及距离半径 $\epsilon$ | | **保守性** | 通常较高，尤其当矩信息范围较宽时 | 可通过 $\theta$ 调节，中等 | 可通过 $\epsilon$ 调节，且具有有限样本概率保证 | | **模型可处理性** | 可转化为半定规划(SDP)或二阶锥规划(SOCP) | 可转化为凸优化问题，有时有解析解 | 常可转化为线性规划(LP)或凸优化，易于求解 | | **优点** | 对数据要求低，仅需矩信息范围 | 直接利用历史数据，概念直观 | 统计性质好，对分布支撑集变化敏感，适合连续分布 | | **缺点** | 忽略高阶矩信息，可能过于保守 | 模糊集仅包含与经验分布支撑相同的分布 | 计算量可能较大，尤其在高维空间 | | **典型应用** | 金融投资组合、设备选址 | 机组组合、能源调度（本文示例） | 最优潮流、库存管理、机器学习 | ```python # 示例：使用Wasserstein距离模糊集的一个简单概念性代码（对比） # 注意：完整实现较复杂，这里仅展示思想。 from scipy.spatial.distance import cdist def wasserstein_distance_1d(p, q, support_p, support_q): """计算一维离散分布p和q之间的1-Wasserstein距离（Earth Mover's Distance）。""" # 假设support已排序 cdf_p = np.cumsum(p) cdf_q = np.cumsum(q) # 通过线性插值，计算在相同分位数点上的逆CDF值，然后求差值的积分 # 简化：对于经验分布，有更高效的计算方法。这里仅为示意。 # 实际中可使用 ot.emd2 来自Python OT库 (POT) pass # 在实际项目中，推荐使用专门的库，如Python的`POT` (Python Optimal Transport)库来计算Wasserstein距离和求解相关DRO问题。 # 安装: pip install pot # 示例代码框架: # import ot # # 计算经验分布 # a = np.ones(N) / N # 历史样本权重 # b = np.ones(M) / M # 可以是另一组样本或均匀分布 # # 计算成本矩阵（样本间的距离） # M_dist = cdist(historical_errors.reshape(-1,1), other_samples.reshape(-1,1)) # # 计算Wasserstein距离 # w_dist = ot.emd2(a, b, M_dist) # 精确计算 # # 基于Wasserstein距离的DRO求解通常涉及求解一个线性规划对偶问题。 ``` 在项目收官阶段，我习惯将完整的数据处理、模型构建和结果分析流程封装成模块化的类或函数，方便复用和参数调优。例如，可以定义一个`WindDROOptimizer`类，其属性包括历史数据、成本参数，方法包括`fit_ambiguity_set`、`solve_single_period`、`solve_multi_period_drcc`等。这样，当新的风电场数据到来时，只需更新数据并调用相应方法即可快速得到鲁棒调度方案。 最后需要强调的是，分布式鲁棒优化不是银弹，它是在**数据有限性**和**决策鲁棒性**之间架起的一座桥梁。它的价值在于提供了一种系统化的方法论，将我们对不确定性的认知（以模糊集的形式）明确地纳入优化模型，从而做出更稳健的决策。在能源转型的浪潮下，随着可再生能源渗透率不断提高，这种能够量化并管理不确定性的工具，其重要性只会日益凸显。