Wasserstein距离的两阶段分布鲁棒优化实践与Matlab实现

markdown复制## 1. 项目概述

### 1.1 核心问题背景

在电力系统调度、供应链管理等实际工程问题中，决策者常常面临双重挑战：既要应对未来可能出现的多种不确定性场景，又要保证决策方案在最坏情况下仍具有可行性。传统随机优化方法依赖精确的概率分布假设，而鲁棒优化又往往过于保守。基于Wasserstein距离的两阶段分布鲁棒优化（DRO）正是为解决这一矛盾而提出的新型建模框架。

我最近在微电网调度项目中实践了这一方法，发现其独特的优势在于：通过Wasserstein球构建的模糊集既能利用历史数据信息，又不需要严格假设数据服从特定分布。下面以Matlab实现为例，详细解析该方法的建模思路与实现技巧。

### 1.2 Wasserstein距离的核心价值

Wasserstein距离（推土机距离）的数学定义为：

$$
W_p(P,Q) = \left( \inf_{\gamma \in \Gamma(P,Q)} \int_{X \times Y} d(x,y)^p d\gamma(x,y) \right)^{1/p}
$$

其中$\Gamma(P,Q)$是边缘分布为$P$和$Q$的所有联合分布集合。在实际应用中，我们通常采用1阶Wasserstein距离（p=1），其具有以下特性：

- **几何直观性**：可解释为将一个分布"搬移"成另一个分布的最小成本
- **数据驱动**：基于有限样本构建的经验分布具有统计一致性
- **计算友好**：对偶问题可转化为线性规划或二次规划

> 注意：Wasserstein半径$\epsilon$的选择至关重要。过小会导致鲁棒性不足，过大则会使解过于保守。实践中建议采用交叉验证法确定。

## 2. 模型构建与对偶转化

### 2.1 两阶段模型框架

典型的两阶段DRO问题可表述为：

$$
\min_{x \in \mathcal{X}} \left\{ c^T x + \sup_{P \in \mathcal{P}} \mathbb{E}_P [Q(x,\xi)] \right\}
$$

其中$Q(x,\xi)$是第二阶段价值函数，$\mathcal{P}$为Wasserstein模糊集：

$$
\mathcal{P} = \{ P : W_1(P,\hat{P}_N) \leq \epsilon \}
$$

### 2.2 对偶转化技术

通过强对偶理论，可将内层的sup问题转化为：

$$
\inf_{\lambda \geq 0} \left\{ \lambda \epsilon + \frac{1}{N} \sum_{i=1}^N \sup_{\xi} [Q(x,\xi) - \lambda \|\xi - \hat{\xi}_i\|] \right\}
$$

这一转化带来三大优势：
1. 将无限维分布问题转为有限维优化
2. 内层sup问题通常可解析求解
3. 外层关于$\lambda$是单变量凸问题

### 2.3 线性决策规则应用

为降低计算复杂度，第二阶段决策常采用仿射策略：

$$
y(\xi) = y_0 + Y \xi
$$

其中$y_0$为基准决策，$Y$为反馈矩阵。这种参数化方式虽然会损失部分最优性，但能显著提升计算效率。

## 3. Matlab实现详解

### 3.1 代码结构设计

```matlab
function [x_opt, obj_val] = Wasserstein_DRO(c, A, b, xi_samples, epsilon)
    % 输入参数：
    % c: 第一阶段成本向量
    % A,b: 第一阶段约束矩阵和向量
    % xi_samples: 不确定性样本矩阵(N×d)
    % epsilon: Wasserstein半径
    
    % 第一阶段决策变量
    x = sdpvar(size(c,1),1);
    
    % 对偶变量
    lambda = sdpvar(1,1);
    
    % 构建对偶问题约束
    constraints = [A*x <= b, lambda >= 0];
    
    % 样本数量
    N = size(xi_samples,1);
    
    % 目标函数构建
    obj = c'*x;
    for i = 1:N
        xi_hat = xi_samples(i,:)';
        [q_val, ~] = solve_second_stage(x, xi_hat, lambda);
        obj = obj + (lambda*epsilon + q_val)/N;
    end
    
    % 求解优化问题
    ops = sdpsettings('solver','gurobi','verbose',0);
    optimize(constraints, obj, ops);
    
    x_opt = value(x);
    obj_val = value(obj);
end

3.2 关键技术实现

3.2.1 第二阶段问题求解

matlab复制function [q_val, y_val] = solve_second_stage(x, xi_hat, lambda)
    % 定义第二阶段决策变量
    y = sdpvar(size(x,1),1);
    
    % 构建目标函数和约束
    q = h(y) - lambda*norm(xi - xi_hat,1);
    constraints = [B*y <= d - C*x, G*xi + H*y <= k];
    
    % 求解优化问题
    optimize(constraints, -q, sdpsettings('verbose',0));
    
    q_val = value(q);
    y_val = value(y);
end

3.2.2 参数调优技巧

1. Wasserstein半径选择：

matlab复制% 交叉验证法选择epsilon
epsilon_grid = linspace(0.1, 2, 10);
cv_loss = zeros(size(epsilon_grid));
for k = 1:length(epsilon_grid)
    cv_loss(k) = cross_validate(train_data, epsilon_grid(k));
end
[~, idx] = min(cv_loss);
epsilon_opt = epsilon_grid(idx);

2. 加速求解策略：

• 使用warm-start初始化
• 对大规模问题采用Benders分解
• 并行化样本循环计算

4. 实战案例：微电网调度

4.1 问题描述

考虑包含光伏、储能和柴油发电机的微电网系统：

• 不确定性：光伏出力预测误差
• 第一阶段决策：日前发电计划
• 第二阶段决策：实时功率调整

4.2 模型参数设置

matlab复制% [成本参数](https://taotoken.net?utm_source=ai)
c_gen = [0.3; 0.5];  % 火电、柴油机发电成本
c_curt = 1.2;        % 光伏弃光惩罚
c_batt = 0.1;        % 储能充放电损耗

% 技术参数
Pmax = [100; 50];    % 机组最大出力
Bmax = 30;           % 储能容量
Ramp = [20; 15];     % 爬坡速率

% 不确定性样本
load pv_error_samples.mat  % 加载历史预测误差数据
epsilon = 0.5;            % Wasserstein半径

4.3 结果分析

通过改变$\epsilon$观察决策变化：

可见随着$\epsilon$增大，系统预留更多备用容量以应对不确定性，体现了鲁棒性的提升。

5. 常见问题与解决方案

5.1 计算效率问题

问题现象：

• 样本数超过1000时求解时间显著增加
• 内存占用过高

解决方案：

1. 采用场景削减技术：

matlab复制% K-means聚类减少样本数
[idx, C] = kmeans(xi_samples, 100);
weights = accumarray(idx,1)/N;

2. 使用稀疏矩阵存储约束
3. 启用Gurobi的Barrier算法并行计算

5.2 保守性控制

问题现象：

• 最优解过于保守导致经济性差

调整策略：

1. 动态调整$\epsilon$：
   $$
   \epsilon(N) = \frac{C}{\sqrt{N}}
   $$
2. 引入风险偏好参数：

matlab复制alpha = 0.9;  % 风险偏好系数
obj = alpha*(c'*x) + (1-alpha)*sup_term;

5.3 数值稳定性

问题现象：

• 对偶问题出现数值震荡
• 收敛困难

改进措施：

1. 添加正则化项：

matlab复制obj = obj + 1e-6*norm(x,2);

2. 约束条件标准化：

matlab复制A_norm = A./max(abs(A),[],2);
b_norm = b./max(abs(A),[],2);

6. 扩展应用方向

在实际项目中，我发现这套方法还可以拓展到以下场景：

1. 供应链库存管理：

• 第一阶段：确定长期采购合约
• 第二阶段：根据实时需求调整短期采购
• 不确定性：市场需求波动

2. 金融投资组合：

• 第一阶段：资产配置比例
• 第二阶段：动态再平衡策略
• 不确定性：资产收益率变化

3. 交通路径规划：

• 第一阶段：出发前路径选择
• 第二阶段：途中实时调整
• 不确定性：路段拥堵情况

对于想要进一步研究的读者，建议从以下文献入手：

• 《Distributionally Robust Optimization: A Review》 (Rahimian & Mehrotra, 2019)
• 《Data-Driven Distributionally Robust Optimization Using the Wasserstein Metric》 (Esfahani & Kuhn, 2018)

实现过程中最关键的体会是：Wasserstein DRO的威力不仅来自数学理论的严谨性，更在于其与实际问题结合的灵活性。通过调整模糊集的形式和决策规则，可以平衡鲁棒性与计算复杂度，这在工程实践中尤为重要。

code复制