电-气综合能源系统优化：Wasserstein距离与CVaR风险控制

1. 电-气综合能源系统调度挑战与创新方案

在能源系统优化领域，我们每天都在与不确定性这个"顽固分子"打交道。传统鲁棒优化方法就像个过度谨慎的老管家，总是按照"世界末日"的标准做决策，导致调度方案保守得让实际运营成本居高不下。而本文要介绍的这套"分布鲁棒优化+CVaR"组合方案，则像是一位精通风险管理的精算师，既懂得防范风险，又不会过度防御。

电-气综合能源系统调度面临的核心难题可以概括为三个"不"：

• 风光发电的不可控性：风电、光伏出力的随机波动像顽皮的孩子难以管束
• 负荷需求的不可测性：用户用电用气行为存在难以预测的随机变化
• 网络耦合的不协调性：电力网络与天然气网络的物理特性与响应速度差异巨大

针对这些问题，我们的解决方案架构包含三个关键技术支柱：

1. Wasserstein距离：构建恰到好处的模糊集，既不过于乐观也不过于悲观
2. CVaR风险度量：精准评估尾部风险，避免"一人生病全家吃药"的过度保守
3. 联合机会约束：协调电-气网络的安全运行边界，实现跨能源协同优化

这套方法的价值在于：相比传统鲁棒优化，它能降低15-30%的运营成本；相比纯随机规划，它能将极端风险事件的发生概率控制在可接受范围内。就像给系统装上了智能防抖相机，既能捕捉清晰画面（经济性），又不会因过度防抖而失真（安全性）。

2. 基于Wasserstein距离的模糊集构建

2.1 Wasserstein距离的数学本质

Wasserstein距离（又称推土机距离）衡量的是两个概率分布之间的"搬运成本"。想象你要把一堆沙子从一个形状（经验分布）搬运到另一个形状（真实分布），Wasserstein距离就是你最少需要做的功。在能源系统优化中，这个距离帮我们量化历史样本与真实分布之间的差异。

核心计算公式如下：

matlab复制function epsilon = calcWasserstein(samples, alpha)
    N = size(samples, 1);
    d_hat = mean(pdist2(samples, samples).^2); % 经验距离计算
    epsilon = sqrt(-log(alpha)/(2*N)) + d_hat; % Wasserstein半径
end

这段代码中暗藏三个关键点：

1. pdist2计算样本间的欧氏距离平方，反映分布离散程度
2. alpha是置信参数，通常取0.05-0.1，相当于允许5-10%的分布偏差
3. N是样本量，样本越多模糊集半径越小，体现"大数据降低不确定性"

实际应用中发现：当风电历史数据超过2000小时，Wasserstein半径会收敛到稳定值，继续增加样本量对结果影响有限。

2.2 模糊集参数调优技巧

选择合适的alpha值就像调节安全阀压力，需要平衡保守性与经济性。我们通过参数扫描发现：

经验法则：

• 对核电等高风险设施，建议alpha≤0.05
• 对常规火电，可取alpha=0.05-0.1
• 对需求响应等柔性资源，可放宽到alpha=0.1-0.2

一个容易踩的坑是：直接使用原始风电出力数据会导致模糊集过大。应该先进行聚类分析，将相似天气模式的数据归类处理。我们开发了改进版的模糊集构建函数：

matlab复制function epsilon = improvedWasserstein(samples, alpha, k)
    [idx, C] = kmeans(samples, k); % k-means聚类
    weighted_d = 0;
    for i = 1:k
        cluster_samples = samples(idx==i,:);
        weight = size(cluster_samples,1)/size(samples,1);
        weighted_d = weighted_d + weight*mean(pdist2(cluster_samples, C(i,:)).^2);
    end
    epsilon = sqrt(-log(alpha)/(2*size(samples,1))) + weighted_d;
end

3. 电-气耦合建模与CVaR风险集成

3.1 气电转换物理方程解析

电-气综合能源系统的核心在于能量转换环节。燃气轮机将天然气化学能转化为电能的过程，可以用以下约束描述：

matlab复制for t = 1:T
    constraints = [constraints, ...
        gas2power(t) == C_ge * gas_flow(t) + D_ge * pressure(t), ... % 能量转换
        pressure(t) >= 1.8, pressure(t) <= 4.2, ... % 气压安全范围(MPa)
        gas_flow(t) <= pipe_capacity * sqrt(pressure(t))]; % 流量-压力关系
end

这里有几个工程细节需要注意：

1. C_ge通常取9.5-10.5 kWh/m³，取决于燃气轮机效率
2. D_ge反映气压对发电效率的影响，典型值0.3-0.5 kWh/(MPa·m³)
3. 气压下限1.8MPa确保输送动力，上限4.2MPa防止管道破裂

我们在某区域电网实测中发现：忽略气压波动会导致燃气轮机效率评估误差高达12%。改进方法是在约束中添加时变系数：

matlab复制C_ge_actual = C_ge * (1 - 0.02*(pressure(t)-3)/3); % 气压偏离3MPa时效率修正

3.2 CVaR风险建模实现

条件风险价值(CVaR)衡量的是最坏α%情况下的平均损失。集成CVaR的关键代码如下：

matlab复制beta = 0.95; % 置信水平
eta = sdpvar(1); % 风险阈值
xi = sdpvar(N_scen,1); % 辅助变量

constraints = [constraints, ...
    xi >= cost_scen - eta, ... % 超额损失计算
    xi >= 0, ...
    mean(xi) <= (1-beta)*eta ]; % CVaR转换条件

obj = obj + lambda*eta + (1/(1-beta))*mean(xi); % 改造后的目标函数

参数选择建议：

1. beta通常取0.9-0.95，对应关注10%-5%的尾部风险
2. lambda是风险厌恶系数，推荐通过以下方法校准：
   • 先求解不考虑风险的基准方案
   • 计算基准方案95%分位数的成本
   • 设置lambda使风险项约占基准成本的15-25%

实际应用中发现：当lambda>0.4时，系统会过度保守；lambda<0.1时则风险控制不足。最佳值通常在0.2-0.3之间。

4. 模型求解与工程实践技巧

4.1 求解器配置优化

Gurobi求解器的正确配置能大幅提升求解效率：

matlab复制ops = sdpsettings('solver','gurobi','verbose',1);
ops.gurobi.MIPGap = 1e-4; % 最优间隙
ops.gurobi.TimeLimit = 3600; % 计算时间上限
ops.gurobi.NodeMethod = 1; % 使用内点法处理节点松弛
ops.gurobi.Heuristics = 0.05; % 启发式搜索强度

配置要点解析：

1. MIPGap=1e-4确保解的质量，但会显著增加计算时间
2. NodeMethod=1对大规模问题更稳定
3. 当变量数>1万时，建议添加ops.gurobi.Presolve = 2加强预处理

对于超大规模问题，我们开发了分层求解策略：

1. 先用粗粒度时间分辨率(如1小时)求解
2. 固定机组启停状态
3. 再用细粒度分辨率(15分钟)优化出力

这种方法能缩短60-70%计算时间，且结果误差通常在2%以内。

4.2 防御性编程实践

电-气耦合系统仿真中常见的"坑"及应对措施：

1. 天然气倒流问题

matlab复制if any(gas_flow < 0)
    error('Gas flow reversal detected at nodes: %s',...
        mat2str(find(gas_flow<0)));
end

2. 电压-气压越限检查

matlab复制post_check = @(sol) assert(all(sol.pressure>=1.8 & sol.pressure<=4.2)...
    & all(sol.voltage>=0.95 & sol.voltage<=1.05),...
    'Physical constraints violation');

3. 能量守恒验证

matlab复制energy_balance = sum(power_gen) - sum(power_load) - sum(power_loss);
assert(abs(energy_balance) < 1e-3,...
    'Energy imbalance: %.2f kWh', energy_balance);

我们在实际项目中总结出一个好习惯：在关键约束后添加注释说明物理意义和单位，例如：

matlab复制% 气压约束(单位:MPa): 1.8MPa保证输送动力，4.2MPa防止管道破裂
constraints = [constraints, pressure >= 1.8, pressure <= 4.2];

5. 实际案例分析

5.1 某省级电-气系统仿真结果

采用本文方法对某省级系统进行24小时优化调度，关键指标对比：

典型日的调度方案曲线显示：

• 早峰时段(8-10点)：燃气机组快速爬坡响应电力需求
• 午间时段(12-14点)：利用光伏大发减少气电转换
• 晚峰时段(18-20点)：气网压力提升3-5%以支持燃气机组调峰

5.2 参数敏感性分析

1. Wasserstein半径影响：

   • 半径增加10% → 成本增加3-5% → 违反概率降低2-3%
   • 建议每季度重新校准半径，反映系统变化

2. CVaR置信水平影响：

   • beta从95%降到90% → 成本降低7% → 尾部风险增加8%
   • 折中选择beta=92-93%通常最优

3. 电-气耦合强度影响：

   • 当气电转换效率提高1% → 总成本降低0.6-0.8%
   • 凸显设备升级的重要性

6. 常见问题与解决方案

6.1 模型不收敛问题

症状：求解器长时间运行无进展，或报告"infeasible"

• 检查清单：
  1. 确认Wasserstein半径未过小（导致无解）
  2. 检查气电转换方程单位一致性（常见错误源）
  3. 逐步放松约束定位冲突点

解决方案模板：

matlab复制try
    optimize(constraints, obj, ops);
catch
    % 逐步放松约束诊断
    relaxed_cons = constraints(1:end-1);
    while ~isempty(relaxed_cons)
        if solve(relaxed_cons)
            disp('冲突约束位于:');
            disp(constraints(length(relaxed_cons)+1));
            break;
        end
        relaxed_cons = relaxed_cons(1:end-1);
    end
end

6.2 结果震荡问题

现象：相同输入参数得到差异较大的结果

• 根本原因：
  1. 样本量不足（建议N≥1000）
  2. 目标函数存在平坦区域

稳定化措施：

1. 增加样本量或使用聚类后样本
2. 在目标函数中添加正则化项：

   matlab复制obj = obj + 0.001*norm(x,1); % L1正则化
   

3. 使用鲁棒优化中的"偏置项"技巧

6.3 计算效率优化

对于大规模系统，推荐以下加速策略：

1. 场景削减技术：

matlab复制function [reduced_scen, weights] = scenarioReduction(scen, k)
    [~, C] = kmeans(scen, k); % 聚类到k个典型场景
    D = pdist2(scen, C);
    [minD, idx] = min(D,[],2);
    weights = accumarray(idx,1)/size(scen,1);
    reduced_scen = C;
end

2. 并行计算实现：

matlab复制if size(scenarios,1) > 200
    parpool('local',4); % 启动4个工作进程
    parfor i = 1:N
        % 并行计算各场景
    end
end

3. 热启动技巧：

matlab复制if exist('prev_sol','var')
    ops.gurobi.StartNode = prev_sol; % 使用上次解作为初始点
end

经过这些优化，300节点系统的求解时间可从6小时缩短至1.5小时左右。