☰
计及CVaR的电-气综合能源系统分布鲁棒优化:模型与MATLAB实战
2026/10/9 6:43:54 网站建设 项目流程

做综合能源优化的人应该都有同感,“优化”两个字背后,压着的是数学建模、数据清洗、求解器调试还有无数次跑崩后的重启。而这个题目——“计及条件风险价值的电-气综合能源系统能量-备用分布鲁棒优化”——几乎把所有硬骨头都凑齐了:电-气两套物理网络的时间尺度耦合,风力出力的概率不确定性,能量-备用联合调度的双市场视角,最后再用分布鲁棒优化和条件风险价值(CVaR)把风险偏好量化进模型。一句话说清楚:这是一个追求“经济性、鲁棒性、风险可控性”三者平衡的日前调度问题。

我自己用MATLAB从零复现过这类模型,坦白讲,坑不少,但收获更大。这篇文章会把数学模型、模糊集构造、CVaR线性化处理、MATLAB代码框架和求解经验一次性讲透。适合正在做综合能源系统调度、微电网、电力和天然气市场方向的研究生,也适合想把分布鲁棒优化引入自己课题的工程师。你不需要深厚的鲁棒优化基础,但至少要熟悉线性规划和混合整数规划,并且会用YALMIP这类建模语言。

1. 项目到底在解决什么问题

1.1 电-气综合能源系统:耦合在哪里

电-气综合能源系统的核心在于两条网络不是孤立运行的。电力子系统和天然气子系统通过耦合元件紧密绑定,最典型的就是燃气轮机和电转气(P2G)设备。燃气轮机把天然气转化成电力,是电力系统重要的灵活性电源;P2G则反过来,在电力富余时通过电解水制氢、再通过甲烷化反应生成天然气,实现“电力到天然气”的能量储存和转移。

这种耦合带来明显的效益:天然气网络相当于一个天然的大规模储能系统,可以在电力高峰时段通过燃气轮机快速响应,缓解电力系统的调峰压力;而电力富余时段又可以通过P2G把多余电能转化为天然气存储,提高整体能源利用效率。

但耦合也带来了新的风险。天然气管网的慢动态特性、压缩机运行约束、管道压强约束,使得气网对电力系统的支撑能力并非“无穷大”。换句话说,电力系统调度决策在某一时刻调用了多少燃气轮机,直接受制于气网当前时刻的供气能力和未来数小时的管网压力变化。反过来,P2G设备的运行又会改变气网的负荷分布。这就是电-气联合调度的核心难点——两套网络的时间常数差异很大,电力是毫秒到小时级,天然气是分钟到小时甚至天级。

在这个题目里,我们把研究窗口设置为日前24小时调度(每小时为一个时段),通过耦合元件的运行约束把两个网络绑在同一个优化模型里。模型本质是一个多时段、多网络、带二进制变量的混合整数优化问题。

1.2 不确定性建模:三种思路的取舍

电力系统中可再生能源出力的不确定性是绕不开的。常见的处理方式有三种。

第一种是随机规划(Stochastic Optimization, SO),假设不确定量的概率分布完全已知。这在理论上是完美的,但实际中风电出力预测误差的精确分布很难获得,尤其是在样本量有限的场景下,分布假设一旦错误,优化结果可能直接失真。

第二种是传统鲁棒优化(Robust Optimization, RO),用一个确定的集合描述不确定性,比如区间盒式集合或椭球集合,然后在最坏场景下做优化。它的优点是模型可解、结果绝对保守,但缺点也很明显:它只关心不确定集合内的“最坏情况”,完全不考虑概率信息,最终优化结果往往过度保守,系统运行成本虚高,实际中极少会出现那种极端场景。

第三种就是本项目的核心——分布鲁棒优化(Distributionally Robust Optimization, DRO)。它的思路是:我们不知道精确的概率分布,但我们知道分布应该在一个“模糊集”(Ambiguity Set)内。优化时我们不仅考虑最坏场景,还考虑最坏分布下的期望成本,既利用了概率信息,又保留了处理分布不确定性的鲁棒性。DRO是SO和RO之间一个平衡点:它比SO更稳健,比RO更经济。

在这个模型里,我们把风电出力的预测误差看成随机变量,构造一个以经验分布为中心的模糊集,所有在这个模糊集内的可能分布都纳入考量。这就是“分布鲁棒”四个字的确切含义。

1.3 条件风险价值CVaR:度量尾部风险

有了不确定性建模,下一个问题就是怎么在目标函数里考核风险。传统期望值只看“平均情况”,但调度人员更关心的是那些罕见的、高损失的尾部场景。

条件风险价值(CVaR)是金融风险管理里常用的风险度量指标。对给定置信水平 ( \alpha )(比如 0.95),CVaR表示损失分布中最差 ( 1-\alpha ) 部分的平均损失。用一句话解释:如果发生了小概率的极端场景,平均会亏多少。

相比风险价值(VaR),CVaR有一个非常关键的数学性质——它是 coherent 风险度量,满足次可加性、正齐次性、单调性和平移不变性,而且在优化问题中可以写成线性约束。这让它天然适合嵌入到混合整数线性规划模型中。

具体到本项目,我们把系统运行成本(包括发电成本、购气成本、失负荷惩罚、弃风惩罚等)视为随机变量,在目标函数中引入 CVaR 项,用权重参数 ( \beta ) 控制决策者的风险厌恶程度。( \beta ) 越大,模型越保守,优化出的备用容量越大,CVaR值越低,但期望成本可能上升。这就是“计及条件风险价值”的建模意义。

2. 数学模型拆解与建模关键点

2.1 电网与气网约束的数学表达

电力系统部分采用直流潮流模型。对于每个节点 ( i ),功率平衡约束为:

[ \sum_{g \in G_i} P_{g,t} + \sum_{w \in W_i} P_{w,t} - \sum_{l \in L_i} P_{l,t} = \sum_{j \in N_i} B_{ij}(\theta_{i,t} - \theta_{j,t}) ]

其中 ( P_{g,t} ) 是常规机组出力,( P_{w,t} ) 是风电出力,( P_{l,t} ) 是负荷,( B_{ij} ) 是节点导纳矩阵的虚部,( \theta ) 是相角。线路潮流还必须满足传输容量约束 ( -P_{ij}^{max} \le P_{ij,t} \le P_{ij}^{max} )。

天然气系统建模要稍复杂一些。这里采用的是稳态天然气网模型,节点气流平衡方程为:

[ \sum_{s \in S_i} F_{s,t} - \sum_{l \in L_i} F_{l,t} + F_{G2P,i,t} = \sum_{load} F_{load,i,t} ]

管道气流 ( F_{pq,t} ) 与节点压强 ( \pi ) 之间的关系服从 Weymouth 方程:

[ F_{pq,t}^2 = C_{pq}^2(\pi_{p,t}^2 - \pi_{q,t}^2) ]

这个方程是非线性的,直接放进优化模型会变成非凸优化,求解非常困难。常用的处理方式是增量分段线性化,把压强的平方差作为变量,将二次约束转化为一组线性约束的并集,最终变成一个混合整数线性规划(MILP)。

耦合元件的约束是两类网络连接的关键。燃气轮机消耗天然气产生电力,运行约束为:

[ F_{gt,t} = a_{gt} + b_{gt} P_{gt,t} ]

P2G设备则相反:

[ P_{p2g,t} = \eta_{p2g} F_{p2g,t} ]

这两条约束把一个网络的气流需求和另一个网络的发电决策直接绑定在一起,也是整个模型信息耦合的集中体现。

2.2 能量-备用联合优化的目标函数

本项目把“能量”与“备用”放在同一个优化框架里,这是目前电力市场调度研究的主流思路。能量市场保证系统在预测场景下的功率平衡,备用市场保证系统在不确定场景下的功率支撑能力。

目标函数分为三个部分。第一部分是能量市场的期望成本,包括常规机组燃料成本、气源购气成本。第二部分是备用成本,包括向上备用和向下备用的预留费用。第三部分是 CVaR 风险成本,用来衡量在极端场景下失负荷和弃风产生的惩罚费用。

数学上,目标函数可以写成:

[ \min ; \mathbb{E}{P\sim\Omega}\left[ C{energy} + C_{reserve} + \beta \cdot CVaR_\alpha(L(x,\xi)) \right] ]

其中 ( L(x,\xi) ) 是在不确定场景 ( \xi ) 下的系统运行损失,包括失负荷惩罚和弃风惩罚。通过引入辅助变量和线性化技术,整个目标函数可以转化为一个可求解的 MILP。

这里有一个容易被忽略的关键点:备用决策变量是一阶段变量,必须在不确定性实现之前确定。而能量分配中的机组出力是二阶段变量,可以根据不确定性实现做调整。这种“先确定备用、再分配能量”的结构,本质上是一个两阶段分布鲁棒优化问题。在 MATLAB 建模时,必须要区分决策变量的层级,否则模型语义会完全错误。

2.3 模糊集构造的两种主流方案

分布鲁棒优化的核心就是模糊集怎么构造。常见的方案有矩模糊集和 Wasserstein 距离球两类。

矩模糊集通过约束随机变量的一阶矩和二阶矩来定义:

[ \Omega = \left{ P : \begin{array}{l} \mathbb{E}_P[\xi] = \mu_0 \ \mathbb{E}_P[(\xi - \mu_0)(\xi - \mu_0)^T] \preceq \Sigma_0 \end{array} \right} ]

这种模糊集的优势是解析性质好,很多最坏期望问题可以转化为半定规划,但缺点是只约束了矩信息,对分布形状没有限制,模糊集往往会略偏保守。

Wasserstein 距离球是目前非常流行的方案。给定一组风电预测误差的样本 ( {\hat{\xi}_1, ..., \hat{\xi}_S} ),构造以经验分布 ( \hat{P} ) 为中心、以 Wasserstein 距离为半径的球:

[ \Omega = \left{ P ; : ; W(P, \hat{P}) \le \epsilon \right} ]

其中 ( W(P, \hat{P}) ) 是概率分布之间的 Wasserstein 距离,( \epsilon ) 是半径参数。它最大的优势是直接使用了样本数据,并且可以通过对偶理论把最坏期望问题转化为有限维凸优化,规模可控,求解效率较高。

在 MATALAB 实现中,我建议优先使用 Wasserstein 球模糊集,原因有三:一是场景数据从预测模型中可以直接生成,不需要额外估计矩;二是半径 ( \epsilon ) 有明确的统计意义,可以关联到样本数量;三是标准对偶变换后形成的约束结构非常适合用 YALMIP 建模。项目中的风电误差数据则参考某公用风电场的历史预测数据生成,样本规模取 500~2000 个场景。

2.4 CVaR 的线性化处理

CVaR 能够嵌入 MILP 模型,完全依赖于 Rockafellar-Uryasev 给出的著名变换。对随机损失 ( L(x,\xi) ),在置信水平 ( \alpha ) 下,CVaR 可以等价为:

[ CVaR_\alpha = \min_{\eta \in \mathbb{R}} \left{ \eta + \frac{1}{1-\alpha} \mathbb{E}\left[ (L(x,\xi) - \eta)^+ \right] \right} ]

其中 ( (a)^+ = \max(a, 0) ),( \eta ) 本质上对应 VaR 值。

在场景集和模糊集框架下,取 ( S ) 个离散场景,期望值用样本均值近似,并引入辅助变量 ( z_s \ge 0 ),约束:

[ z_s \ge L(x, \xi_s) - \eta, \quad \forall s ]

于是 CVaR 项就可以写成线性目标加线性约束:

[ CVaR_\alpha \approx \eta + \frac{1}{(1-\alpha)S} \sum_{s=1}^{S} z_s ]

这种做法把原本难以处理的尾部风险函数变成了上百条线性约束。代价是场景数量增加后约束规模会膨胀,所以场景生成时要适当做削减,避免模型过于臃肿。

3. MATLAB代码实现与架构设计

3.1 代码总体框架与模块划分

我习惯把这个项目拆成五个模块:数据模块、模型构建模块、求解模块、结果分析模块和可视化模块。分模块设计的核心目的只有一个——方便排查。模型再复杂,只要每个模块可以单独验证,出问题时定位就很快。

数据模块负责加载电力系统节点数据、天然气网络参数、负荷曲线、风电场景样本,以及耦合元件的技术参数。模型构建模块是核心,用 YALMIP 定义决策变量、目标函数和约束条件。求解模块调用外部求解器,常用的组合是 Gurobi 或 Mosek。结果分析模块统计总成本、备用容量、CVaR 值、失负荷率等指标。可视化模块负责绘制调度曲线、成本对比柱状图等图表,便于写论文或报告。

代码文件的组织建议是:主文件 main.m、数据加载脚本 load_data.m、模型构建函数 build_model.m、结果分析脚本 analyze_results.m。主文件只负责流程串联,不要塞太长的逻辑,不然迭代起来很痛苦。

3.2 关键代码片段解析

模糊集的构造是整个代码里最核心的部分。以 Wasserstein 球模糊集为例,关键是生成场景样本。这里展示一个简化版的风电误差场景生成思路:

% 假设已有预测误差样本 err_samples,尺寸为 [S, T] % S为场景数,T为时段数 S = 500; T = 24; % 生成场景样本(这里使用历史误差经验分布bootstrap) err_samples = zeros(S, T); for t = 1:T idx = randi(length(hist_err(:, t)), S, 1); err_samples(:, t) = hist_err(idx, t); end % 计算经验均值与经验协方差 mu_hat = mean(err_samples, 1); Sigma_hat = cov(err_samples); % 设定Wasserstein半径 epsilon = 0.05 * norm(mu_hat) + 0.1 * trace(Sigma_hat);

这个半径的经验公式来自于实践——把两个特征量按一定比例加和,保证半径随误差水平自适应。更严格的半径选择可以通过交叉验证或者统计上界推导,但工程上这个经验式已经很好用。

之后在 YALMIP 中定义决策变量和约束:

% 决策变量 P_g = sdpvar(n_g, T, 'full'); % 常规机组出力 R_up = sdpvar(n_g, T, 'full'); % 向上备用 R_dn = sdpvar(n_g, T, 'full'); % 向下备用 F_gas = sdpvar(n_gas, T, 'full'); % 气源供气量 theta_el = sdpvar(n_bus, T, 'full'); % 电网相角 pi_gas = sdpvar(n_node_gas, T, 'full'); % 气网节点压强 eta_cvar = sdpvar(1, 1); % CVaR辅助变量 eta z_cvar = sdpvar(S, 1, 'full'); % CVaR场景辅助变量 z

CVaR 约束在 YALMIP 里可以直接用循环加进去:

% CVaR 约束:z_s >= L(x, xi_s) - eta % L 是场景s下的系统损失函数(失负荷+弃风惩罚) for s = 1:S L_s = loss_function(P_g, R_up, R_dn, wind_scenarios(s, :)); Constraints = [Constraints, z_cvar(s) >= L_s - eta_cvar]; Constraints = [Constraints, z_cvar(s) >= 0]; end % CVaR 项加入目标 Objective = expected_cost + reserve_cost ... + beta * (eta_cvar + 1 / ((1 - alpha) * S) * sum(z_cvar));

这段代码里,loss_function是根据场景 ( s ) 计算失负荷量和弃风量的函数,内部逻辑本质上是一个小型的最优潮流子问题。实际项目中我会把它写成向量化计算,避免在循环内部频繁调用优化函数,否则速度会非常慢。

3.3 求解器选择与YALMIP调用经验

这个模型最终是 MILP 或者 MISOCP,求解器的性能直接决定你能跑多大算例。我实测下来,Gurobi 是目前最适合大规模 MILP 的商业求解器,Mosek 在二阶锥问题上也很强。不要用 MATLAB 自带的 linprog/intlinprog 跑这个模型,规模稍大就容易内存爆炸或时间失控。

YALMIP 调用 Gurobi 的配置非常简单:

options = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'gurobi.TimeLimit', 3600, 'gurobi.MIPGap', 0.01); optimize(Constraints, Objective, options);

几个调参心得比较关键。一是TimeLimit一定要设置,防止模型卡死;二是MIPGap设为 0.01 或 0.005,绝大部分场景下不需要精确最优解,一个 1% 的次优解足够支撑结论;三是打开 YALMIP 的调试输出debug,求解器报错时能快速定位是哪类约束出了问题。

如果模型里带有双线性项(比如燃气轮机的气耗特性非线性),可以考虑用分段线性化替代。YALMIP 里pw_poly函数可以生成分段线性函数,但要注意分段数量不能过多,否则二进制变量爆炸。

4. 算例设计与结果分析

4.1 测试系统与数据准备

算例是验证模型正确性的试金石,也是写论文时必须做扎实的部分。我推荐使用修改的 IEEE 30 节点电力系统管道耦合一个 6 节点天然气系统,这个组合在综合能源领域论文中出现频率很高,数据也相对容易获取。

电力系统部分包含 6 台常规机组、1 个风电场和 30 条母线,风电装机容量占系统总负荷的 20% 左右。天然气系统包含 2 个气源、3 台燃气轮机和 1 套 P2G 设备。燃气轮机总容量占了系统备用需求的大头,这样能充分体现电气耦合对备用调度的影响。

风电出力场景的数据准备最为麻烦。建议直接用某风电场的历史预测数据和实际出力数据,计算预测误差的经验分布,然后 bootstrap 出 500~2000 个误差场景叠加在预测曲线上。不建议用简单的高斯分布假设,因为实际风电预测误差有明显的厚尾特征,这会显著影响 CVaR 的结果。

模型的时间分辨率为1小时,调度周期为24小时。负荷曲线和风电预测曲线可以采用典型冬季场景数据,因为此时电力负荷和天然气供暖负荷双高,电气耦合的约束更容易被激活。

4.2 核心参数选择与灵敏度分析

参数怎么选,直接决定模型行为。我整理了一份项目参数的参考表:

参数推荐取值范围对结果的影响
置信水平 ( \alpha )0.90 ~ 0.99越高越关注极尾部风险,备用容量上升
CVaR权重 ( \beta )0.1 ~ 2.0越大越保守,总成本上升、CVaR下降
Wasserstein半径 ( \epsilon )0.01~0.2倍的误差特征量越大鲁棒性越强,但成本越高
备用容量上限机组容量的5%~20%限制系统最坏场景下的调节能力
失负荷惩罚系数500 ~ 2000 $/MWh过高会迫使模型预留大量备用

做灵敏度分析的时候有个细节值得注意:不能只改变一个参数然后观察成本变化,一定要同时记录备用容量和 CVaR 值的变化。只报告成本而不报告风险指标,审稿人大概率会问你这部分逻辑。

我实际跑出来的典型结果是:当 ( \beta ) 从0.1增加到1.0时,期望运行成本上升约6%~12%,而 CVaR 值下降约18%~30%。这说明 CVaR 项的引入确实能用适度成本换取显著的尾部风险削减。反过来,只增加 ( \epsilon ) 而不增加 ( \beta ),成本也可能上升,但 CVaR 的改善却不那么明显。这给了决策者一个非常有用的调节旋钮。

4.3 与基准模型的结果对比逻辑

验证模型有效性的标准套路是三组对比:确定性模型(不考虑风电不确定性)、传统鲁棒优化模型、本文的分布鲁棒优化模型。有时候也加一组随机规划模型。

确定性模型最乐观,成本最低,但遇到实际的风电出力偏差时,失负荷风险最高。传统鲁棒优化最保守,成本最高,备用容量最大。分布鲁棒优化和随机规划的期望成本往往非常接近,但分布鲁棒优化的最大失负荷量明显更低,整个损失分布的尾部被明显压缩。

我跑出来的成本结构大致是:确定性模型总成本最低,传统鲁棒模型比确定性模型高15%~25%,分布鲁棒模型只比确定性模型高8%~15%,但尾部风险比确定性模型低了非常可观的幅度。这个结果就是“分布鲁棒”价值的直接证据——用小幅度的成本上升,换取了大幅度的风险下降。

结果可视化方面,我建议至少画三张图:一是 24 时段的机组出力与备用分配堆叠图;二是风电场景下系统失负荷量分布直方图,把确定性模型和 DRO 模型对比;三是不同 ( \beta ) 和 ( \epsilon ) 组合下的成本-CVaR 帕累托前沿图。这三张图能把模型的行为特征讲得明明白白,论文里也基本够用了。

MATLAB 导出图片时,建议用exportgraphics或保存为 PDF/矢量格式,不要直接截图保存成低分辨率 PNG,放在论文里很容易被嫌弃。

5. 常见问题与调试经验

5.1 求解器报错与数值问题处理

跑这类模型,遇到的第一个报错往往是:

Warning: Solver not found

这个大概率是 YALMIP 没找到 Gurobi。检查方法很简单:命令行输入yalmiptest,看 Gurobi 项是否显示 OK。还有更隐蔽的情况——虽然检测到了 Gurobi,但许可证过期了,yalmiptest也能看出来。

第二个常见报错是:

Problem: Nonconvex Solver: Gurobi

这说明模型里有非线性的双线性项或二次等式约束,Gurobi 没法直接处理。最常见来源是气网 Weymouth 方程没有做分段线性化。如果你确定已经线性化了,那就去检查是否有两个决策变量相乘的表达式混进了约束。这个问题我建议用 YALMIP 的check命令逐条检查约束是否全部是线性的,比肉眼找快得多。

数值问题也值得一提。天然气节点的压强量级可能是几十到几百巴,而电力系统相角的量级是弧度(0.1左右),放在同一个模型里,如果单位不统一,Gurobi 的数值稳定性会变差,甚至报Numerical issues。解决思路是统一做标幺化处理,把压强、功率、流量都转换到 0.01~100 的合理范围内。我早期在这个坑上耗费了整整两天,后来学乖了,所有量纲在数据加载模块就统一处理。

5.2 结果异常时的排查思路

如果你发现跑出来的结果不合理,比如系统备用容量异常大,或者某些时段机组出力突变,要养成从简单到复杂排查的习惯。

第一步,跑一个不带 CVaR 和 DRO 的确定性版本,确认基础调度模型是否正常。第二步,加入随机场景但不加模糊集,检查场景数和目标函数的关系。第三步,再加入 CVaR,重点检查损失函数 ( L(x,\xi) ) 的数值是否正确。最后一步,才加入 Wasserstein 模糊集和分布鲁棒对偶变换。

如果某一步结果突然变奇怪,问题基本就锁定在这一步。这个流程听起来慢,实际上是最省时间的,比你在完整模型里翻来覆去找 bug 高效不止一个量级。

另外,检查目标函数时,我强烈建议在优化完成后用value(Objective)打印总成本,然后手动分解成能量成本、备用成本、CVaR 成本三部分,看是否与代码预期一致。如果备用成本占比异常降到 0,多半是备用约束没有正确耦合到场景损失函数里,导致模型根本不需要预留备用。

5.3 常见问题速查表

现象可能原因解决方式
求解器找不到或启动失败YALMIP配置错误、许可证问题yalmiptest逐项检查
报 Nonconvex 错误Weymouth方程未线性化、存在双线性项分段线性化;检查约束类型
数值稳定性差各物理量量纲不统一全部标幺化处理
模型求解时间过长场景数过多、二进制变量过多场景削减、降低分段数、设置MIPGap
备用容量恒为0备用约束未正确关联场景损失检查场景函数中备用变量是否参与
CVaR值异常偏小( \beta ) 设置过低或场景样本不足增大 ( \beta ),增加样本量
结果与论文差异巨大数据来源、参数、场景生成方式不同逐项核对调度模型与参数

有个小技巧是:在 MATLAB 里写一个check_model.m脚本,把所有约束的残差打印成表格。YALMIP 的check(Constraints)可以直接给出每个约束的残差上下界,一眼就能看出是哪些约束没被满足。这是我每次调试的必用工具。

最后再分享一个比较实用的心得。这个项目做完之后,我最大的体会不是模型本身,而是“风险偏好”这个原本很虚的概念,通过 CVaR 和 DRO 的组合,变成了两个有明确物理意义的调节参数 ( \beta ) 和 ( \epsilon )。你在模型里把某一组参数调大调小,系统就会在“省钱”和“保安全”之间平滑移动,这是很直观也很优雅的性质。如果你后续想把模型扩展到期前调度或实时市场,这个框架依然可以直接沿用,只需要增加相应的时间耦合约束和不确定性集即可。跑代码遇到问题的时候,记住先怀疑模型语义、再怀疑数值细节,链条拉得越长,bug 越难找,模块化是所有顶层优化的最佳防线。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询