Python实战:用Wasserstein距离搞定风光发电预测的分布鲁棒优化(附完整代码)

如果你正在新能源领域,尤其是风光发电预测和电力系统调度中摸爬滚打,那么“不确定性”这个词一定让你又爱又恨。风力和光伏出力天生具有间歇性和波动性,依赖历史数据训练的传统预测模型,在面对极端天气或数据分布偏移时,往往表现不佳,导致调度计划失准,甚至影响电网安全。我们需要的,是一种能坦然承认“未来无法精确预知”,但依然能做出“最坏情况下也不至于太差”的决策方法。这正是分布鲁棒优化(Distributionally Robust Optimization, DRO)大显身手的舞台。

而要让DRO从理论公式落地为工程代码,关键在于如何量化“不确定性”。我们不再假设风光出力服从某个完美的正态分布,而是承认,真实分布可能在我们观测到的历史数据(经验分布)周围“摇摆”。如何划定这个“摇摆”的范围?这就需要一把尺子——Wasserstein距离。它不像KL散度那样对分布的形态有苛刻要求,而是直观地衡量了“将一个分布搬运成另一个分布所需的最小工作量”,非常适合刻画基于有限样本时我们对真实分布认知的模糊地带。

本文将彻底抛开复杂的理论推导,聚焦于一个目标:手把手带你用Python,从零构建一个基于Wasserstein距离的分布鲁棒优化模型,并将其应用于风光发电预测后的调度决策问题。我们将从数据标准化、Wasserstein球半径计算,到最终的鲁棒优化模型构建与求解,提供完整的、可运行的代码块,并穿插可视化和实战技巧,让你不仅能看懂,更能直接用起来。

1. 环境准备与核心概念落地

在开始写代码之前,我们需要确保工具箱是齐全的。除了经典的数值计算库,凸优化求解器是本次实战的核心。

pip install numpy pandas matplotlib scipy cvxpy

这里重点介绍一下cvxpy。它是一个用于构建和求解凸优化问题的Python库,语法非常直观。在DRO问题中,经过对偶转化后,我们常常会得到一个凸优化问题,cvxpy能让我们像写数学公式一样自然地构建模型。

Wasserstein距离与模糊集:工程视角的理解 抛开严格的数学定义,你可以这样理解:我们有一组历史风光出力观测数据,构成了一个“经验分布”。我们承认真实分布可能不是它,但相信真实分布不会离它“太远”。Wasserstein距离就是衡量“远近”的那把尺子。我们设定一个容忍半径ε,所有与经验分布Wasserstein距离不超过ε的分布,就构成了我们的“模糊集”(Ambiguity Set)。DRO的目标,就是在这个模糊集里最坏的那个分布下,寻找最优决策。这保证了只要真实分布落在我们划定的圈子里,我们的决策就不会崩盘。

注意:Wasserstein距离的计算本身可能很复杂,但在构建模糊集时,一个关键步骤是确定这个半径ε。我们将采用基于样本的统计方法进行估计,后文会给出具体代码。

2. 数据预处理与Wasserstein球半径计算

任何数据驱动的模型都始于数据。假设我们有一个风光电站的历史出力数据集samples,其形状为(n_features, m_samples)。例如,n_features=2可以分别代表风电和光伏的归一化出力,m_samples是历史时刻数。

2.1 数据标准化与中心化

标准化不是为了美观,而是为了后续计算Wasserstein距离时,不同量纲的特征具有可比性。我们采用常见的“中心化+缩放”处理。

import numpy as np

def standardize_samples(samples):
    """
    对风光出力样本进行标准化处理。
    参数:
        samples: numpy数组,形状 (n_features, m_samples)
    返回:
        thet: 标准化后的样本矩阵
        mu: 各特征的经验均值
        sig: 各特征的经验标准差
    """
    n, m = samples.shape
    # 计算均值和标准差
    mu = np.mean(samples, axis=1, keepdims=True)  # 形状 (n, 1)
    sig = np.std(samples, axis=1, keepdims=True)   # 形状 (n, 1)
    # 避免除零,给标准差一个极小值
    sig = np.where(sig < 1e-10, 1.0, sig)
    # 标准化: (样本 - 均值) / 标准差
    thet = (samples - mu) / sig
    return thet, mu, sig

# 示例数据:假设有1000个历史时刻,风电和光伏两种出力
np.random.seed(42)
m_samples = 1000
wind_power = np.random.weibull(2.0, m_samples) * 100  # 模拟风电
pv_power = np.random.beta(2, 5, m_samples) * 50       # 模拟光伏
samples_raw = np.vstack([wind_power, pv_power])

thet, mu, sig = standardize_samples(samples_raw)
print(f"原始样本形状: {samples_raw.shape}")
print(f"标准化后样本形状: {thet.shape}")
print(f"风电均值: {mu[0,0]:.2f}, 标准差: {sig[0,0]:.2f}")
print(f"光伏均值: {mu[1,0]:.2f}, 标准差: {sig[1,0]:.2f}")

2.2 计算Wasserstein球半径ε

这是连接经验数据与模糊集的关键桥梁。半径ε控制着模糊集的大小:ε越大,考虑的分布不确定性越强,结果越保守(鲁棒);ε越小,则越接近传统的随机优化(基于经验分布)。一种常见的方法是使用统计学习理论中的边界,其公式通常涉及样本量、置信水平和一个与数据分布相关的常数C。

下面的代码实现了一种基于样本极差(已知支撑集)或经验矩生成函数(未知支撑集)的半径估计方法。我们以实现后者为例,它更通用。

from scipy.optimize import minimize

def compute_wasserstein_radius(thet, beta=0.95, support_known=False, support_range=None):
    """
    计算Wasserstein模糊集的半径 epsilon。
    参数:
        thet: 标准化后的样本矩阵,形状 (n, m)
        beta: 置信水平,通常取0.9, 0.95等
        support_known: 是否已知随机变量的支撑集(取值范围)
        support_range: 如果已知,提供每个特征的取值范围,形状 (n, 2)
    返回:
        epsilon: 计算得到的Wasserstein球半径
    """
    n, m = thet.shape
    # 计算样本的2-范数距离矩阵(这里简化计算,使用样本与均值差的范数)
    # 更严谨的做法是计算所有样本对之间的欧氏距离,但计算量大。这里采用一种近似。
    # 计算每个样本到原点(因为thet已中心化,原点即均值点)的欧氏距离
    dist_to_center = np.linalg.norm(thet, axis=0)  # 形状 (m,)
    
    if support_known and support_range is not None:
        # 如果已知支撑集,半径与支撑集直径相关
        # 假设支撑集是超立方体,计算其直径
        diameters = support_range[:, 1] - support_range[:, 0]  # 每个特征的取值范围长度
        D = np.sqrt(np.sum(diameters**2))  # 超立方体空间对角线长度
        C = D / 2.0  # 一种常用的常数估计
    else:
        # 未知支撑集,通过优化一个目标函数来估计常数C
        # 目标函数: J(alpha) = sqrt( (1/(2*alpha)) * (1 + log( (1/m) * sum(exp(alpha * dist**2) ) ) )
        def obj_c(alpha):
            if alpha <= 0:
                return 1e10  # 加一个很大的惩罚,确保alpha为正
            test = dist_to_center ** 2
            inner_sum = np.mean(np.exp(alpha * test))
            J = np.sqrt( (1/(2*alpha)) * (1 + np.log(inner_sum)) )
            return J
        
        # 使用优化器寻找使J最小的alpha
        res = minimize(obj_c, x0=1.0, method='L-BFGS-B', bounds=[(1e-6, None)])
        alpha_opt = res.x[0]
        C = 2 * obj_c(alpha_opt)
    
    # 计算最终半径 epsilon
    epsilon = C * np.sqrt( (2/m) * np.log(1/(1-beta)) )
    return epsilon

# 计算示例
epsilon = compute_wasserstein_radius(thet, beta=0.95, support_known=False)
print(f"计算得到的Wasserstein球半径 epsilon: {epsilon:.4f}")

这个epsilon值就是我们的“安全边际”。它告诉我们,基于现有的1000个样本,在95%的置信水平下,我们认为真实分布不会离经验分布(标准化后)超过epsilon的Wasserstein距离。

3. 构建风光预测的分布鲁棒优化模型

现在进入核心部分。假设我们是一个微电网的调度员,需要在风电、光伏出力不确定的情况下,决定从主网购电功率P_grid(决策变量),以最小化总成本(购电成本 + 惩罚成本),同时满足负荷需求。

问题设定:

  • 已知一个固定的负荷需求 P_load
  • 风光联合出力 P_renewable 是一个二维随机变量(风电,光伏),其真实分布未知,但属于以经验分布为中心、Wasserstein半径为epsilon的模糊集。
  • 我们可以从主网以单价 c_grid 购电。
  • 如果总供电(购电+风光出力)小于负荷,会产生缺电惩罚,单价为 c_penalty(远高于购电单价)。
  • 目标:在最坏的可能分布下,最小化期望总成本。

这是一个典型的两阶段分布鲁棒优化问题。第一阶段决定购电量(“此时此地”的决策),第二阶段在不确定性揭示后承担惩罚(“看到再看”的补偿)。经过理论推导(利用Wasserstein DRO的对偶变换),该问题可以转化为一个确定性的凸优化问题,具体来说是一个二阶锥规划(SOCP)问题。

import cvxpy as cp

def build_and_solve_dro_scheduling(thet, mu, sig, epsilon, P_load, c_grid, c_penalty):
    """
    构建并求解基于Wasserstein DRO的风光调度问题。
    参数:
        thet, mu, sig: 标准化后的样本及其统计量
        epsilon: Wasserstein球半径
        P_load: 负荷需求 (标量)
        c_grid: 单位购电成本
        c_penalty: 单位缺电惩罚成本
    返回:
        P_grid_opt: 最优购电功率
        objective_value: 最优目标函数值(最坏情况下的期望成本上界)
        status: 求解状态
    """
    n, m = thet.shape  # n: 风光特征数, m: 样本数
    
    # 第一阶段决策变量:从主网的购电功率
    P_grid = cp.Variable(nonneg=True)
    
    # 辅助变量,用于刻画最坏情况分布下的期望惩罚成本
    # 根据对偶理论,我们需要引入变量 lambda, s_i
    lambda_var = cp.Variable(nonneg=True)  # 对偶变量,与epsilon相关
    s = cp.Variable(m)                     # 每个样本对应的辅助变量
    
    # 构建目标函数:购电成本 + 最坏情况下的期望惩罚成本上界
    # 目标函数形式: c_grid * P_grid + lambda * epsilon + (1/m) * sum(s_i)
    objective = cp.Minimize(c_grid * P_grid + lambda_var * epsilon + (1/m) * cp.sum(s))
    
    # 构建约束条件
    constraints = []
    
    # 对于每一个样本 i,需要满足的约束(源自对偶变换)
    for i in range(m):
        # 该样本对应的风光出力(还原到原始尺度)
        # xi = mu + sig * thet[:, i],这里thet[:, i]是标准化后的扰动
        xi = mu + sig * thet[:, i].reshape(-1, 1)  # 形状 (n, 1)
        # 总可再生能源出力(假设风电和光伏出力可直接相加)
        P_renew_i = cp.sum(xi)  # 标量
        
        # 缺电量(如果购电+可再生小于负荷)
        # 注意:在优化模型中,我们直接定义缺电惩罚项为 hinge loss: max(0, P_load - P_grid - P_renew_i)
        # 根据对偶理论,这导出一个关于s_i和lambda的约束
        # 约束形式: s_i >= c_penalty * [P_load - P_grid - P_renew_i] - lambda * ||thet[:, i]||_2
        # 并且 s_i >= 0
        shortage = P_load - P_grid - P_renew_i
        norm_theta_i = cp.norm(thet[:, i], 2)  # 样本i的2-范数
        constraints.append(s[i] >= c_penalty * shortage - lambda_var * norm_theta_i)
        constraints.append(s[i] >= 0)
    
    # 购电功率上限约束(可选)
    P_grid_max = P_load * 1.5  # 假设最大购电能力为负荷的1.5倍
    constraints.append(P_grid <= P_grid_max)
    
    # 定义问题并求解
    prob = cp.Problem(objective, constraints)
    # 使用ECOS或SCS求解器,它们擅长处理二阶锥约束(cp.norm隐含了)
    try:
        prob.solve(solver=cp.ECOS, verbose=False)
    except Exception as e:
        print(f"ECOS求解失败: {e},尝试使用SCS求解器。")
        prob.solve(solver=cp.SCS, verbose=False)
    
    # 获取结果
    if prob.status in ['optimal', 'optimal_inaccurate']:
        P_grid_opt = P_grid.value
        obj_val = prob.value
        status = prob.status
        print(f"求解成功!状态: {status}")
        print(f"最优购电功率: {P_grid_opt:.2f} kW")
        print(f"最坏情况下期望总成本上界: {obj_val:.2f} 元")
        
        # 计算基于经验分布的成本(作为对比)
        # 即假设未来风光出力就是历史样本的等概率混合,计算平均成本
        P_renew_samples = np.sum(mu + sig * thet, axis=0)  # 每个样本的总可再生出力
        cost_samples = c_grid * P_grid_opt + c_penalty * np.maximum(0, P_load - P_grid_opt - P_renew_samples)
        empirical_cost = np.mean(cost_samples)
        print(f"基于经验分布的平均成本: {empirical_cost:.2f} 元")
        print(f"鲁棒优化付出的‘保守溢价’: {obj_val - empirical_cost:.2f} 元")
        
        return P_grid_opt, obj_val, status
    else:
        print(f"求解失败,状态: {prob.status}")
        return None, None, prob.status

# 设置问题参数
P_load = 200.0  # kW
c_grid = 1.0    # 元/kWh
c_penalty = 10.0 # 元/kWh (惩罚远高于购电成本)

# 调用函数求解
P_opt, obj_val, status = build_and_solve_dro_scheduling(thet, mu, sig, epsilon, P_load, c_grid, c_penalty)

这段代码构建了整个DRO调度模型的核心。关键在于理解对偶变换后引入的变量lambda_varslambda_var实质上权衡了不确定性的大小(ε)与惩罚成本,而s_i捕获了每个样本场景下的成本。约束条件 s_i >= c_penalty * shortage - lambda * ||theta_i|| 正是Wasserstein DRO理论中对偶变换的具体体现,它将一个无限维的分布优化问题,转化为了一个有限维的凸优化问题。

4. 模型分析与可视化:鲁棒性如何体现?

代码跑通了,结果也出来了。但这个“鲁棒”的调度方案,到底比传统方法好在哪里?我们需要更直观的感受。

4.1 对比不同优化策略

我们来对比三种策略:

  1. 确定性优化:假设风光出力等于其历史均值,直接求解。
  2. 随机优化(基于经验分布):假设未来风光出力的分布就是历史样本的等概率混合,最小化期望成本。
  3. 分布鲁棒优化(本文方法):考虑以经验分布为中心、Wasserstein球内的最坏分布。

我们将计算在不同真实分布(模拟一些可能偏离历史经验的情况)下,三种策略的实际成本。

import matplotlib.pyplot as plt

def evaluate_policy(P_grid_decided, true_samples, P_load, c_grid, c_penalty):
    """评估给定购电决策在实际分布下的平均成本。"""
    P_renew_true = np.sum(true_samples, axis=0)
    costs = c_grid * P_grid_decided + c_penalty * np.maximum(0, P_load - P_grid_decided - P_renew_true)
    return np.mean(costs)

# 1. 确定性优化(按均值决策)
P_renew_mean = np.sum(mu)  # 风光出力的均值
P_grid_det = max(0, P_load - P_renew_mean)  # 只需补充缺额
print(f"\n确定性优化决策: 购电 {P_grid_det:.2f} kW")

# 2. 随机优化(基于经验分布)—— 等价于上述DRO模型中令 epsilon=0
# 我们可以通过调用DRO函数,但设置epsilon为一个极小的值来近似
epsilon_tiny = 1e-6
P_grid_sto, obj_sto, _ = build_and_solve_dro_scheduling(thet, mu, sig, epsilon_tiny, P_load, c_grid, c_penalty)

# 3. 分布鲁棒优化(使用之前计算的epsilon)
P_grid_dro = P_opt

# 生成几种可能的“真实分布”进行测试
np.random.seed(123)
test_costs = {'确定性': [], '随机优化': [], '分布鲁棒': []}
n_test_scenarios = 5

for i in range(n_test_scenarios):
    # 模拟真实分布偏离经验分布:例如,风光同时偏弱(均值下移)
    shift_factor = np.random.uniform(-0.3, 0.1, size=(2, 1))  # 风电和光伏分别的偏移系数
    scale_factor = np.random.uniform(0.8, 1.2, size=(2, 1))   # 波动性变化
    
    # 生成测试样本
    true_mu = mu * (1 + shift_factor)
    true_sig = sig * scale_factor
    # 假设测试分布也是正态的(仅为演示,实际可能非正态)
    m_test = 5000
    true_thet = np.random.randn(2, m_test)  # 标准正态
    true_samples = true_mu + true_sig * true_thet
    
    # 评估三种策略
    cost_det = evaluate_policy(P_grid_det, true_samples, P_load, c_grid, c_penalty)
    cost_sto = evaluate_policy(P_grid_sto, true_samples, P_load, c_grid, c_penalty)
    cost_dro = evaluate_policy(P_grid_dro, true_samples, P_load, c_grid, c_penalty)
    
    test_costs['确定性'].append(cost_det)
    test_costs['随机优化'].append(cost_sto)
    test_costs['分布鲁棒'].append(cost_dro)

# 可视化对比
fig, ax = plt.subplots(figsize=(10, 6))
x = np.arange(n_test_scenarios)
width = 0.25
ax.bar(x - width, test_costs['确定性'], width, label='确定性优化', color='skyblue')
ax.bar(x, test_costs['随机优化'], width, label='随机优化(经验分布)', color='lightgreen')
ax.bar(x + width, test_costs['分布鲁棒'], width, label='分布鲁棒优化', color='salmon')

ax.set_xlabel('测试场景编号')
ax.set_ylabel('平均实际成本 (元)')
ax.set_title('不同优化策略在不同真实分布下的性能对比')
ax.set_xticks(x)
ax.legend()
ax.grid(True, axis='y', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()

# 计算平均成本和最坏情况成本
print("\n=== 跨场景性能总结 ===")
for policy in test_costs:
    costs = test_costs[policy]
    print(f"{policy}:")
    print(f"  平均成本: {np.mean(costs):.2f} 元")
    print(f"  最坏场景成本: {np.max(costs):.2f} 元")
    print(f"  成本标准差: {np.std(costs):.2f} 元")

4.2 敏感性分析:Wasserstein半径ε的影响

半径ε是控制模型保守程度的核心“旋钮”。ε=0时,模型退化为传统的随机优化(完全信任经验分布);ε越大,模型越保守,考虑的不确定性越强。我们来观察最优购电决策和最优成本如何随ε变化。

# 分析不同epsilon对决策和成本的影响
epsilon_list = np.linspace(0.0, epsilon*2.0, 15)  # 从0到2倍的计算半径
P_grid_list = []
obj_value_list = []

for eps in epsilon_list:
    P_opt_eps, obj_val_eps, status_eps = build_and_solve_dro_scheduling(thet, mu, sig, eps, P_load, c_grid, c_penalty)
    if P_opt_eps is not None:
        P_grid_list.append(P_opt_eps)
        obj_value_list.append(obj_val_eps)
    else:
        P_grid_list.append(np.nan)
        obj_value_list.append(np.nan)

# 绘制结果
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(epsilon_list, P_grid_list, 'o-', linewidth=2, markersize=8)
ax1.set_xlabel('Wasserstein球半径 ε')
ax1.set_ylabel('最优购电功率 (kW)')
ax1.set_title('最优决策随ε的变化')
ax1.grid(True, linestyle='--', alpha=0.7)
ax1.axvline(x=epsilon, color='r', linestyle='--', label=f'计算半径 (ε={epsilon:.3f})')
ax1.legend()

ax2.plot(epsilon_list, obj_value_list, 's-', linewidth=2, markersize=8, color='orange')
ax2.set_xlabel('Wasserstein球半径 ε')
ax2.set_ylabel('最坏情况期望成本上界 (元)')
ax2.set_title('最优成本上界随ε的变化')
ax2.grid(True, linestyle='--', alpha=0.7)
ax2.axvline(x=epsilon, color='r', linestyle='--', label=f'计算半径 (ε={epsilon:.3f})')
ax2.legend()

plt.tight_layout()
plt.show()

通过这张图,你可以清晰地看到模型的“保守性-经济性”权衡。随着ε增大,系统为了防范更恶劣的风光出力场景,会倾向于购买更多的主网电力作为备用,导致购电成本上升(决策更保守),同时最坏情况成本上界也随之提高。在实际应用中,ε的选取可以基于历史数据的丰富程度和对风险的厌恶程度来调整。

5. 进阶讨论与工程实践要点

将理论模型投入实际生产环境,还需要考虑更多细节。以下是几个关键点:

1. 模糊集的构建不止Wasserstein距离一种 Wasserstein模糊集因其良好的统计性质和易处理性而流行,但它不是唯一选择。其他常见类型包括:

模糊集类型 核心思想 优点 缺点
矩模糊集 约束真实分布的前几阶矩(如均值、协方差)在一定范围内。 模型通常可转化为半定规划(SDP),有成熟求解器。 仅利用矩信息,可能过于保守,且高阶矩估计不准。
φ-散度模糊集 使用KL散度、χ²散度等衡量分布差异。 对某些问题可导出线性或锥优化问题。 要求分布绝对连续,对支撑集敏感。
Wasserstein模糊集 基于分布之间的“搬运”距离。 具有良好的数据驱动特性,有限样本性能有理论保证。 计算可能更复杂,对偶问题规模与样本数相关。

选择哪种模糊集,取决于具体问题的特性、数据的质量以及可接受的计算复杂度。

2. 处理高维与大规模样本 我们的示例中风光特征只有2维。在实际电网中,可能需要考虑数十个甚至上百个风电场/光伏电站,此时thet的维度n会很大。此外,历史样本数m也可能达到数万级别。这会导致优化问题中约束数量急剧增加(我们代码中的循环约束有m个)。

应对策略

  • 场景削减:使用K-means等聚类方法,将m个历史样本聚合成K个典型场景(K << m),用聚类中心代表一类样本,可以大幅减少约束数量。
  • 分布式计算:对于大规模问题,可以考虑使用分布式优化算法,或者利用cvxpy结合高性能求解器(如MOSEK)的并行特性。
  • 采样近似:使用随机梯度下降(SGD)或随机近似方法求解DRO的对偶问题,适用于样本数极大的情况。

3. 与机器学习预测模型的结合 本文假设风光出力的历史样本是直接可用的。但在实际中,我们通常使用机器学习模型(如LSTM、Transformer)进行点预测或概率预测。DRO可以很自然地与这些模型结合:

  • 点预测+误差分布:用模型预测出力值,将历史预测误差作为thet的样本。DRO用于刻画预测误差分布的不确定性。
  • 概率预测输出:如果模型能输出预测分布(如分位数、参数化分布),可以将该预测分布作为DRO的“中心分布”,并围绕其构建模糊集。

4. 代码优化与生产部署建议

  • 模型封装:将DRO模型构建、求解、结果后处理封装成一个类,方便参数管理和调用。
  • 求解器选择:对于中小规模问题,ECOSSCS是不错的选择。对于更大规模或需要更高精度的问题,可以考虑商业求解器如MOSEKGurobi
  • 缓存与热启动:如果只是参数(如epsilon, c_penalty)微调,而问题结构不变,可以利用上一次求解的结果作为本次求解的初始点,加速收敛。
  • 结果验证:始终在历史数据或留出的测试集上回测策略,并与简单的基准策略(如确定性规则)比较,确保鲁棒性提升不是以牺牲过多经济性为代价。
# 一个简单的模型封装示例
class WassersteinDROScheduler:
    def __init__(self, historical_samples, P_load, c_grid, c_penalty, beta=0.95):
        self.samples_raw = historical_samples
        self.P_load = P_load
        self.c_grid = c_grid
        self.c_penalty = c_penalty
        self.beta = beta
        self.thet, self.mu, self.sig = None, None, None
        self.epsilon = None
        self._fitted = False
        
    def fit(self):
        """拟合模型:标准化数据并计算Wasserstein半径。"""
        self.thet, self.mu, self.sig = standardize_samples(self.samples_raw)
        self.epsilon = compute_wasserstein_radius(self.thet, beta=self.beta)
        self._fitted = True
        print(f"模型拟合完成。计算得到 epsilon = {self.epsilon:.4f}")
        
    def solve(self, epsilon=None):
        """求解DRO调度问题。"""
        if not self._fitted:
            raise ValueError("请先调用 .fit() 方法拟合模型。")
        if epsilon is None:
            epsilon = self.epsilon
            
        P_opt, obj_val, status = build_and_solve_dro_scheduling(
            self.thet, self.mu, self.sig, epsilon, self.P_load, self.c_grid, self.c_penalty
        )
        self.solution_ = {'P_grid': P_opt, 'objective': obj_val, 'status': status}
        return self.solution_
    
    def evaluate(self, test_samples):
        """在测试样本集上评估当前决策的性能。"""
        if self.solution_ is None:
            raise ValueError("请先调用 .solve() 方法获得决策。")
        P_grid = self.solution_['P_grid']
        avg_cost = evaluate_policy(P_grid, test_samples, self.P_load, self.c_grid, self.c_penalty)
        return avg_cost

# 使用示例
scheduler = WassersteinDROScheduler(samples_raw, P_load=200, c_grid=1.0, c_penalty=10.0, beta=0.95)
scheduler.fit()
solution = scheduler.solve()
print(f"调度决策: 购电 {solution['P_grid']:.2f} kW")

风光预测的不确定性管理是一个持续的过程。分布鲁棒优化提供了一套严谨的数学框架,让我们能在承认认知局限的前提下,做出更具韧性的决策。本文提供的代码和思路是一个起点,你可以根据具体的应用场景调整模糊集的形式、目标函数和约束。例如,在电力系统中,你可能还需要考虑线路潮流约束、储能设备动作约束等,这些都可以整合到cvxpy的建模框架中。关键在于理解Wasserstein距离如何将数据的不确定性转化为优化问题中一个可调节的“保守度”参数,从而在过度乐观和极端保守之间找到属于你当前系统的最佳平衡点。

更多推荐