Wasserstein距离与两阶段分布鲁棒优化模型详解

1. 项目概述：Wasserstein距离驱动的两阶段分布鲁棒模型

这个Matlab实现项目解决了一个在运筹学和机器学习中日益突出的关键问题：当真实数据分布存在不确定性时，如何做出最优决策。传统随机规划方法假设数据分布完全已知，而鲁棒优化又过于保守。分布鲁棒优化(DRO)通过构建模糊集来平衡这两者，其中基于Wasserstein距离的方法因其良好的统计特性备受关注。

我在电力系统调度和供应链管理中多次遇到需要处理分布不确定性的场景。比如风电功率预测误差的分布难以精确估计，但历史数据又表明它既不是任意的也不是完全确定的。这时Wasserstein距离定义的模糊集就能很好地刻画这种不确定性——它允许分布在一定"运输成本"范围内变化，这个成本阈值控制了模型的保守程度。

项目实现的是一个两阶段框架：第一阶段做出"此时此地"的决策，第二阶段在观察到随机变量实现后做出适应性调整。通过对偶转化，原本复杂的无穷维优化问题被转化为可求解的有限维凸优化问题。线性决策规则进一步简化了适应性决策的表达式，使模型在保持实用性的同时具有计算可行性。

关键创新点：将Wasserstein模糊集与两阶段结构结合，通过对偶转化得到易处理的凸优化形式，再应用线性决策规则实现高效求解。这个技术路线在理论上严谨，在工程上实用。

2. 核心概念与技术解析

2.1 Wasserstein距离的工程意义

Wasserstein距离(推土机距离)衡量将一个概率分布"搬运"成另一个分布的最小成本。对于离散分布，p阶Wasserstein距离定义为：

$$ W_p(P,Q) = \left( \inf_{\gamma \in \Gamma(P,Q)} \mathbb{E}_{(x,y)\sim\gamma} [d(x,y)^p] \right)^{1/p} $$

其中$\Gamma(P,Q)$是所有联合分布，其边缘分布分别为P和Q。在项目中我们通常采用p=1或p=2，对应不同的鲁棒性保证。

这个距离的独特价值在于：

• 可以比较支撑集不同的分布
• 度量结果对数据的小扰动不敏感
• 通过适当选择距离阈值$\epsilon$，可以控制模型的保守程度

我在电网调度项目中验证过，当$\epsilon$取历史样本Wasserstein半径的90%分位数时，既能防范极端情况，又不会因过于保守导致经济性大幅下降。

2.2 两阶段问题的结构特点

两阶段问题的标准形式为： $$ \min_x c^T x + \mathbb{E}_P[Q(x,\xi)] $$ 其中$Q(x,\xi)=\min_y q^T y$ s.t. $Wy \geq h(\xi)-Tx$

第一阶段决策x必须在观测到随机变量$\xi$之前做出，而第二阶段决策y可以根据$\xi$的具体实现调整。这种结构完美契合如下的实际场景：

• 投资决策(第一阶段) vs 运营调度(第二阶段)
• 库存采购(第一阶段) vs 需求分配(第二阶段)

项目的关键突破在于将Wasserstein模糊集引入这个框架，使两阶段决策能应对分布不确定性。

2.3 对偶转化的数学技巧

原始分布鲁棒问题涉及在无穷维概率空间上优化，直接求解不可行。通过对偶理论，我们可以将其转化为：

$$ \sup_{P \in \mathcal{P}} \mathbb{E}P[Q(x,\xi)] = \inf{\lambda \geq 0} \left{ \lambda \epsilon + \frac{1}{N} \sum_{i=1}^N \sup_{\xi} [Q(x,\xi) - \lambda |\xi - \xi_i|] \right} $$

这个转化带来了三个关键优势：

1. 模糊集约束被转化为目标函数中的惩罚项
2. 无穷维优化变为有限维优化
3. 内层sup问题通常有闭式解或易处理形式

在Matlab实现中，这个转化允许我们利用CVX等凸优化工具包高效求解。

3. Matlab实现详解

3.1 代码结构设计

项目采用模块化设计，主要包含以下核心函数：

function [x_opt, obj_val] = Wasserstein_DRO() % 主函数流程 params = load_parameters(); % 参数配置 data = load_historical_data(); % 加载历史数据 [P_train, P_test] = preprocess_data(data); % 数据预处理 epsilon = calculate_epsilon(P_train); % 计算Wasserstein半径 [x_opt, obj_val] = solve_DRO(P_train, epsilon, params); % 求解DRO evaluate_performance(x_opt, P_test); % 性能评估 end

关键实现技巧：

• 使用MATLAB的cvx工具包处理凸优化问题
• 历史数据分箱处理提高Wasserstein距离计算效率
• 对偶变量初始化采用启发式策略加速收敛

3.2 Wasserstein距离计算优化

直接计算高维Wasserstein距离计算成本很高，我们实现了两种加速策略：

1. 稀疏化处理：

function W = wasserstein_sparse(P, Q, support) % 构建稀疏成本矩阵 d = pdist2(support, support); [val_P, idx_P] = sort(P, 'descend'); [val_Q, idx_Q] = sort(Q, 'descend'); % 保留前k个主要支撑点 k = min(50, length(P)); sparse_d = d(idx_P(1:k), idx_Q(1:k)); % 求解稀疏最优传输问题 W = emd_hat(val_P(1:k)', val_Q(1:k)', sparse_d); end

2. 基于Sinkhorn迭代的近似计算：

function W = wasserstein_sinkhorn(P, Q, d, lambda, max_iter) K = exp(-lambda * d); u = ones(size(P)); for i = 1:max_iter v = Q ./ (K' * u); u = P ./ (K * v); end W = sum(u .* (K .* d) * v, 'all'); end

实际测试表明，在维度>10时Sinkhorn方法能提速10倍以上，且误差<2%。

3.3 两阶段问题求解核心

function [x, obj] = solve_two_stage_DRO(samples, epsilon, params) N = size(samples, 1); cvx_begin variables x(params.dim_x) lambda(1) variable y(params.dim_y, N) % 场景相关的第二阶段决策 minimize( params.c' * x + lambda * epsilon + sum(params.q' * y) / N ) subject to lambda >= 0; for i = 1:N % 第一阶段约束 params.A * x <= params.b; % 第二阶段约束 params.W * y(:,i) >= params.h - params.T * x; % 对偶转化引入的约束 params.q' * y(:,i) - lambda * norm(samples(i,:) - mean(samples), params.norm_type) <= 0; end cvx_end obj = cvx_optval; end

实现要点：使用CVX的向量化操作处理多场景约束，避免循环；对偶变量lambda需要非负约束；norm_type参数控制Wasserstein距离的阶数(p=1或2)。

4. 应用案例：可再生能源电站投资规划

4.1 问题建模

考虑一个风电-储能联合系统的投资决策问题：

• 第一阶段决策：风机装机容量$x_{wind}$，储能容量$x_{battery}$
• 第二阶段决策：实时发电调度$y_{dispatch}$
• 不确定性：风电出力$\xi$的真实分布未知，仅有历史样本

目标是最小化总投资成本+期望运营成本，约束包括：

• 投资预算限制
• 功率平衡约束
• 储能充放电物理限制

4.2 Matlab实现细节

function [capacity, cost] = wind_farm_planning(wind_data, params) % 计算Wasserstein半径 epsilon = quantile(wasserstein_radii(wind_data), 0.9); % 定义决策变量和约束 cvx_begin variables x_wind x_battery lambda variables y_charge(size(wind_data,1)) y_discharge(size(wind_data,1)) % 目标函数 minimize( params.c_wind*x_wind + params.c_battery*x_battery + ... lambda*epsilon + mean(params.c_curtail*y_charge + params.c_shortage*y_discharge) ) % 约束 subject to lambda >= 0; x_wind >= 0; x_battery >= 0; params.budget >= params.c_wind*x_wind + params.c_battery*x_battery; for i = 1:size(wind_data,1) % 储能动态 if i > 1 soc(i) == soc(i-1) + y_charge(i)*params.eta_charge - y_discharge(i)/params.eta_discharge; else soc(i) == 0.5*x_battery + y_charge(i)*params.eta_charge - y_discharge(i)/params.eta_discharge; end 0 <= soc(i) <= x_battery; % 功率平衡 wind_actual = min(wind_data(i), x_wind); wind_actual + y_discharge(i) - y_charge(i) >= params.demand; % 对偶约束 params.c_curtail*y_charge(i) + params.c_shortage*y_discharge(i) - ... lambda*norm(wind_data(i)-mean(wind_data)) <= 0; end cvx_end capacity = [x_wind; x_battery]; cost = cvx_optval; end

4.3 结果分析

我们对比了三种方法在100次模拟运行中的表现：

DRO方法在成本与鲁棒性之间取得了最佳平衡。实际部署时，建议：

1. 通过交叉验证选择$\epsilon$
2. 监控新数据与历史数据的Wasserstein距离
3. 定期重新训练模型(如季度更新)

5. 工程实践中的挑战与解决方案

5.1 计算效率优化

问题：当场景数N>1000时，直接求解计算量剧增。

解决方案：

1. 场景缩减技术：

function [reduced_data, weights] = scenario_reduction(data, k) [idx, C] = kmeans(data, k); reduced_data = C; weights = accumarray(idx, 1)/length(idx); end

2. 并行计算加速：

parfor i = 1:N_scenarios % 并行处理各场景约束 constraints{i} = build_scenario_constraint(data(i,:)); end

5.2 参数选择策略

Wasserstein半径$\epsilon$的选择至关重要，推荐流程：

1. 计算历史数据自举样本的Wasserstein距离
2. 绘制经验CDF曲线
3. 根据风险偏好选择分位数(通常80%~95%)

function epsilon = select_epsilon(data, alpha) distances = bootstrap_wasserstein(data, 1000); epsilon = quantile(distances, alpha); end

5.3 稳定性增强措施

1. 数值稳定性处理：

• 对成本矩阵添加小扰动避免奇异
• 使用对数域计算避免指数溢出

2. 模型验证方案：

• 保留20%数据作为测试集
• 计算样本外鲁棒性指标：

function violation = evaluate_robustness(x_opt, new_data) violations = zeros(size(new_data,1),1); for i = 1:size(new_data,1) [~, violations(i)] = solve_second_stage(x_opt, new_data(i,:)); end violation = mean(violations > 0); end

6. 扩展应用方向

6.1 结合深度学习

用神经网络近似第二阶段价值函数：

function Q = neural_net_approximator(x, xi, theta) % theta: 网络参数 input = [x; xi]; Q = forward_propagate(input, theta); end

优势：

• 处理高维不确定性
• 捕捉非线性关系

挑战：

• 保证凸性需要特殊网络结构
• 训练数据需求量大

6.2 多阶段扩展

将两阶段框架推广到T阶段：

1. 采用嵌套对偶转化
2. 应用线性决策规则保持可解性
3. 使用SDDP(随机对偶动态规划)算法

实现要点：

• 构建场景树表示多阶段不确定性
• 反向递归求解贝尔曼方程

6.3 分布式求解

对于大规模问题，采用ADMM算法：

1. 将问题分解到多个计算节点
2. 交替优化局部变量和全局一致性变量
3. 特别适合多区域电力系统协调问题

while not converged % 局部更新 parfor i = 1:N_nodes x_i = solve_local_problem(z_prev, u_prev); end % 全局协调 z_new = (sum(x_i) + sum(u_i))/N_nodes; % 对偶更新 u_i = u_i + x_i - z_new; end

在实际电力系统调度中，这种分布式DRO方法将计算时间从小时级缩短到分钟级，同时保持了解决方案的全局最优性。