1. 项目背景为什么偏偏是“电-气综合能源系统”和“能量-备用”一起优化我做综合能源系统优化也有几年了最开始接触的其实只有电力系统那时候觉得经济调度无非就是让机组出力最小化约束就那几条。后来参与了一个园区级电-气耦合项目才意识到单独看电网根本不够——天然气管道、压缩机、储气罐和燃气机组耦合在一起任何一个环节的波动都会传导到电力侧。这次分享的“计及条件风险价值的电-气综合能源系统能量-备用分布鲁棒优化”本质上就是在解决一个问题当新能源出力和负荷预测都存在误差时怎么同时安排电能量出力、气源供气以及向上/向下备用才能既保证经济性又控制极端场景下的风险。这里的“能量”包括电功率平衡、天然气流量平衡、气源产量、燃气机组耗气量这些“备用”指的是系统应对预测偏差而预留的调节能力比如燃煤机组上调出力、可中断负荷、储气罐释放等。传统做法通常是能量调度和备用调度分开做或者备用只按固定比例预留但电-气耦合系统里备用响应会直接影响天然气流量需要联动优化。项目标题里的“条件风险价值”和“分布鲁棒优化”是两个撬动核心问题的支点。CVaR用来度量“最差的一部分场景”的平均损失分布鲁棒则允许我们在只知道部分统计信息如一阶矩、二阶矩的情况下构造一个包含真实分布的不确定集然后求解最坏情况下的期望成本。两者结合既不像传统随机优化那样依赖精准分布假设也不像传统鲁棒优化那样只盯着最坏场景而过于保守。这个思路放在电-气综合能源系统里特别适合处理新能源出力预测误差的非参数不确定性和备用需求的随机性。对于刚接触这个方向的同学我的建议是先把“为什么不能简单用期望值”想清楚。用风电举例某时刻预测出力100MW但实际可能只有60MW或者140MW。如果调度时只按100MW安排电量一旦实际偏低就要大量调用高价备用如果按最坏情况60MW安排又会因为少用便宜风电而浪费经济性。分布鲁棒优化就是在这两者之间找一个平衡通过模糊集控制“最坏情况下的平均损失”本质上是给决策者的风险偏好加上一个可调参数。2. 分布鲁棒优化的建模思想为什么要用Wasserstein距离构造模糊集2.1 从不确定性建模说起随机与鲁棒的折中先梳理一下常见的不确定性处理方式。随机优化假设随机变量的概率分布完全已知比如正态分布然后求解期望成本最小化。问题很简单但真实的风电出力并不服从标准正态尤其是尾端常常比正态更厚用正态分布可能严重低估极端情况。传统鲁棒优化则假设不确定量落在一个确定的集合内比如盒式集合或预算不确定集目标是在最坏情况下满足约束结果往往异常保守运行成本居高不下。分布鲁棒优化DRO是两者的中间路线我们不知道精确分布但知道一些统计特征或者有一组历史样本于是构造一个“候选分布集合”模糊集真实分布以大概率落在这个集合内然后求解最坏情况下的期望成本。这个“最坏情况下的期望”跟CVaR天然有联系因为CVaR本身就是关于某个分布下尾部损失的条件期望。2.2 模糊集的选择对比为什么我倾向Wasserstein球模糊集有多种构造方式常见的有矩约束模糊集、Kullback-Leibler散度球、Wasserstein球等。矩约束模糊集直接用均值和方差的波动范围来约束建模简单但对分布形状刻画弱。KL散度球需要基准分布且只适用于连续分布间距离度量数值求解时对支撑集有要求。Wasserstein距离则基于最优运输思想可以度量两个分布之间“物质搬移”的最小成本因为在几何意义上更直观而且对分布形状没有太多限制近年来在综合能源优化里很受欢迎。我用的是Wasserstein球模糊集以经验分布为中心半径为ε的Wasserstein球内所有分布都被视为可能候选。半径ε越大模糊集越大决策越保守ε0时就退化为样本均值近似即随机优化的经验版本。在实际程序中ε通常通过历史样本数量、置信度和调度时长来确定但做项目时也可以把它当作风险偏好参数来调节。2.3 CVaR与分布鲁棒的结合点条件风险价值CVaR定义为在置信水平α下损失超过VaR值的尾部损失的平均值。用公式表达就是给定随机变量L比如运行成本CVaR_α(L) E[L | L ≥ VaR_α(L)]。在优化模型里一般引入辅助变量η和松弛变量来线性化变成CVaR_α(L) η (1/(1-α)) E[ max(L - η, 0) ]这个形式在MATLAB中用YALMIP的cvar命令或自己写辅助变量都能实现。将CVaR嵌入分布鲁棒目标函数时我们不再是单纯最小化最坏情况期望而是最小化“最坏情况下超过某个门槛的平均损失”。这样既可以控制高风险场景又避免了传统鲁棒那种“一刀切”式的极端优化。在电-气综合能源系统里备用容量不足导致的切负荷损失、弃风损失、气源中断损失都可以统一折算成货币成本而这些损失往往集中在极端天气或故障场景中用CVaR来测度非常合适。3. 电-气综合能源系统的能量-备用联合建模3.1 电力系统子模型要点电力侧我用的是直流潮流模型忽略无功和电压只关心有功平衡。节点功率平衡约束为∑_i P_{G,i,t} ∑_w P_{W,w,t} - ∑_j P_{L,j,t} ∑_b B_{bus}·θ_t这里P_G是常规机组出力P_W是新能源出力比如风电P_L是负荷B_bus是节点导纳矩阵的虚部B矩阵θ是节点相角。直流模型对系统级调度够用误差在可接受范围内计算效率高。备用约束包括上备用和下备用∑_{i∈G} r_{G,i,t}^{up} ∑_{s∈S} r_{S,s,t}^{up} ≥ R_t^{up}简单理解就是所有能上调出力的资源机组、储能、可中断负荷、气转电之和要满足系统上调备用需求。下调备用同理。备用需求R_t^{up}本身不是固定值而是和新能源预测误差、负荷预测误差相关联这正是分布鲁棒优化的切入对象。3.2 天然气系统子模型要点天然气网络通常用Weymouth方程描述管道流量与节点压力的非线性关系。为了便于求解很多文献用分段线性化或二阶锥松弛。我这次用的是分段线性化因为我和YALMIP里的binvar配合比较熟好控制精度。气源节点有产量上限储气罐有充放气速率和容量约束压缩机消耗部分天然气但同时能提升下游压力。针对电-气耦合的核心是燃气机组耗气量F_{g,t} a·P_{g,t}^2 b·P_{g,t} c这个二次关系我同样做了分段线性化。另一个关键约束是燃气机组的备用容量会转化为额外耗气需求也就是说上备用意味着燃气机组可能会多发电需要气网提供额外流量。这是电-气耦合备用优化的最主要体现——传统独立调度往往忽略这条导致事故发生时气源供气不足。3.3 能量与备用的协同决策变量设计决策变量分两层第一层日前调度常规机组出力、气源产气量、储气罐计划、备用容量预留。第二层实时再调度当随机实现发生后各资源调整出力但调整量必须被限制在日前预留的备用容量内。这样本质上是一个两阶段分布鲁棒优化min_{x∈X} c·x sup_{P∈D} E_P[ Q(x, ξ) ]其中Q(x, ξ)是第二阶段最优成本包括高成本备用调用、切负荷、弃风、气源短缺等惩罚。用CVaR时第二阶段成本的目标被替换为CVaR度量。这个结构在MATLAB里可以直接用嵌套优化不行需要用对偶或分解算法。我采用的方式是将第二阶段线性规划对偶后合并到主问题把分布鲁棒部分通过对偶转化转化为有限维凸优化。这个转化细节比较繁琐我在下面章节详细说。4. MATLAB实现从数学公式到可运行代码4.1 编程环境与关键工具我用的版本是MATLAB R2024a优化求解器是Gurobi 10通过YALMIP作为建模层。没有Gurobi的话用Mosek或者CPLEX也行只要有二次锥规划能力。因为我用的分段线性化和Wasserstein球对偶后都是线性和二阶锥约束Gurobi非常合适。注意不要直接使用MATLAB自带的linprog因为模型规模一大求解速度会明显下降。典型环境配置% 添加YALMIP路径然后调用Gurobi求解器 addpath(genpath(C:\yalmip)); ops sdpsettings(solver,gurobi,verbose,2,gurobi.MIPGap,0.01);4.2 数据初始化与系统参数以IEEE 30节点电力系统和6节点天然气系统耦合为例。电力系统部分有6台常规机组、2个风电场、3台燃气机组天然气系统部分有2个气源、3条管道、1个储气罐。基础参数初始化% 基础参数 T 24; % 调度时段 nb 30; % 电力节点数 ng 6; % 常规机组数 ngas 3; % 燃气机组数 nwind 2; % 风电场数 nload 20; % 负荷节点数 np 6; % 天然气节点数 nsource 2; % 气源数 npipe 3; % 管道数 % 负荷预测曲线 load_pred(i,t)风电预测 wind_pred(w,t) % 历史误差样本 hist_wind_error(s,t)用于构造模糊集 % 常规机组成本系数 c_cost, c_reserve_up, c_reserve_dn % 切负荷惩罚系数 c_curtail, c_shed 等历史误差样本通常从历史运行数据中提取比如过去30天每天24个时段的风电预测误差和负荷预测误差。如果没有真实数据可以用正态分布加厚尾模拟生成但代码注释里要标明这是模拟数据。4.3 两阶段分布鲁棒建模核心代码建模的关键在于第二阶段的决策变量和约束怎么表达。我直接展示核心结构% 第一阶段变量 P_G sdpvar(ng, T, full); % 常规机组出力 R_up sdpvar(ng, T, full); % 上备用预留 R_dn sdpvar(ng, T, full); % 下备用预留 Q_source sdpvar(nsource, T, full); % 气源产气 V_store sdpvar(nstore, T, full); % 储气罐状态 % 第二阶段变量针对每个场景 P_G_fix sdpvar(ng, T, full); % 调整后的出力 % 实际场景中风电出力 预测值 误差 % 调整量受限于备用 P_G_fix(i,t) - P_G(i,t) R_up(i,t); P_G_fix(i,t) - P_G(i,t) -R_dn(i,t);如果直接把场景数设为几千变量会爆炸。所以我不展开场景而是用分布鲁棒对偶把第二阶段期望部分转化为线性矩阵不等式。具体来说对于Wasserstein球模糊集第二阶段目标可以写成min_{y∈Y(x,ξ)} d·y其最坏期望问题可以用强对偶转化为sup_{P∈D} E_P[ Q(x,ξ) ] inf_{λ≥0} λ·ε (1/N)∑_{k1}^N sup_{ξ∈Ξ} [ Q(x,ξ) - λ·c(ξ, ξ_k) ]其中c(ξ, ξ_k)是Wasserstein距离的代价函数。当支撑集Ξ是有限维凸集时内部sup问题可以进一步线性化。最终可以在YALMIP里写成lambda sdpvar(1); % Wasserstein半径对偶变量 radius 0.5; % 可调模糊集半径 objective c_first*[P_G; R_up; R_dn; Q_source] lambda*radius (1/N)*sum(eta_k); % eta_k 第二阶段最优值 - lambda*距离函数具体公式推导在下一篇帖子里详细展开MATLAB里推荐直接用cvar的YALMIP扩展来简化CVaR写法但底层还是需要自己建模。4.4 目标函数与CVaR的YALMIP表示目标函数由三部分组成第一阶段成本机组发电成本、气源产气成本、备用预留成本、储气成本。第二阶段CVaR由于不确定性导致的再调度成本、切负荷惩罚、弃风惩罚。如果天然气网无法满足燃气机组额外用气还会加气源短缺惩罚。YALMIP代码可以这样组织% 第一阶段成本 cost_first sum(C_G(:).*P_G(:)) sum(C_Rup(:).*R_up(:)) sum(C_Rdn(:).*R_dn(:)); cost_source sum(C_source(:).*Q_source(:)); cost_store sum(C_store(:).*abs(V_store(:)-V_store0(:))); % 第二阶段损失函数构造 loss_fun (xi) second_stage_cost(xi, x_decision); % x_decision是第一阶段变量xi是随机向量风电误差、负荷误差 % 分布鲁棒CVaR目标 obj cost_first cost_source cost_store distributionally_robust_cvar(loss_fun, samples, alpha, epsilon_w); ops sdpsettings(solver,gurobi); optimize(cons, obj, ops);这里的distributionally_robust_cvar函数需要自己实现。最直接的方式是引入辅助变量eta_k和t_k把每一段距离惩罚项线性化最终交给Gurobi。重点代码eta sdpvar(1); % CVaR阈值 z_k sdpvar(N,1); % 尾端损失辅助变量 % 对每个样本k约束 z_k 0 且 z_k loss_k - eta % loss_k是对应样本的最坏再调度成本由内层max给出 % 最终CVaR eta (1/(1-alpha))* (1/N)*sum(z_k)4.5 如何将非线性约束线性化电-气模型主要有三个非线性来源第一管道流量Weymouth方程F_pipe k * sqrt(p_up^2 - p_down^2)我用分段线性化逼近。做法是把流量范围切成M段每段对应一组二进制变量逐段约束精度控制在2%以内。对于系统级调度误差可以接受。第二燃气机组耗气量F_g a P_g^2 b P_g c用分段线性插值。用interp1生成分段系数。第三储气罐的库容状态方程包含压差项但我采用简化的线性动态模型V(t1)V(t)Q_in(t)-Q_out(t)这在园区级规划尺度上是合理的。线性化之后的模型是混合整数线性规划加二阶锥Gurobi处理起来很快。5. 算例设置与结果分析从对比里看清模型价值5.1 测试系统参数我测试的耦合系统是修改后的IEEE 30节点比利时20节点天然气网络。不过为了突出电-气互动我简化成6节点气网参数如下参数数值风电预测误差标准差预测值的15%负荷预测误差标准差预测值的5%置信水平α0.95Wasserstein半径ε0.5 (归一化)备用成本倍率发电成本的30%切负荷惩罚5000 $/MWh弃风惩罚200 $/MWh模糊集所需历史样本设为100组样本从误差历史序列中随机抽样生成。5.2 与随机优化、传统鲁棒的对比我分别实现三个对照组样本平均近似SAA、传统盒式鲁棒优化、本文的分布鲁棒CVaR方法。算例结果如下成本单位$10^5方法预期运行成本最坏场景成本成本标准差计算时间SAA12.318.72.445s传统鲁棒16.814.91.160s本文DROCVaR13.515.81.6110s注意最坏场景成本这一列传统鲁棒的最低因为它在设计时已经强制要求任何情况下都必须满足约束但代价是预期运行成本高得离谱。SAA虽然预期成本最低但最坏场景成本高达18.7说明它在面对未知分布时非常脆弱。本文方法预期成本比SAA高9.7%但最坏场景成本降低了15.5%标准差也小得多。这个结果说明分布鲁棒CVaR的折中特性非常明显牺牲一点平均经济性换来尾部风险显著下降。5.3 备用配置结果的关键洞察看调度结果时我特别关注了一个典型日风电预测大误差时段14:00-16:00。SAA给出的备用容量很小因为它假设预测误差服从正态没有为厚尾场景留余地。传统鲁棒给出的备用量几乎让所有机组都处于热备用状态显然不经济。本文方法则把备用主要配置在两部分一是位于电-气耦合负载率较高区域的燃气机组上二是储气罐释放能力上。原因很简单燃气机组上备用时天然气管道可能会成为瓶颈。如果只在电力侧保留备用当燃气机组实际需要追加出力时气网供气跟不上备用形同虚设。分布鲁棒优化通过CVaR识别了“天然气供应不足”这一尾部风险所以会在气源侧额外预留在产气能力同时给储气罐设置较高的初始充放气裕度。这种跨系统备用配置是单一能量系统优化做不到的。5.4 不同Wasserstein半径和置信水平的影响我还做了敏感性分析。半径ε从0.2到2.0变化运行成本从12.8升到15.2最坏场景成本从17.1降到14.7。这个趋势跟理论一致半径越大模糊集越大决策越保守。置信水平α从0.9到0.99变化时CVaR对尾部损失越敏感备用预留越多。实际操作中建议根据运行人员对风险的忍受程度来选参数没有绝对最佳。置信水平α影响的是“尾部”的定义α0.95意味着只关心最差的5%场景如果系统安全标准要求苛刻可以调到0.99。而Wasserstein半径α则可以理解为对环境不确定性的信任程度样本量大、数据质量高时可以选小值数据稀疏时选大值。6. 常见坑与调试心得亲手踩过的雷全部分享6.1 数值爆炸问题我第一次跑电-气耦合的分布鲁棒优化直接用了单层optimize()结果Gurobi报错说“Numerical trouble”。排查后发现是分段线性化系数差异太大——电力侧成本系数在10^2量级切负荷惩罚在10^3量级而气网压力单位换算后是10^6量级。这些系数量级差距过大导致KKT矩阵病态。解决思路就是归一化。我把所有成本统一折算成标幺值以100MVA为基准压力用标幺值。比如管道压力取基准值6.3MPa流量基准取200kmol/h。归一化之后求解器再没报过数值错误。这个细节如果你用我的代码数据一定要提前处理。6.2 YALMIP中CVaR辅助变量的约束错误很多新手在写CVaR时会把eta也当作变量加上下限实际上eta在标准CVaR公式中是无界变量但优化问题必须有界才可解。所以需要给eta加一个大范围边界比如[-1e6, 1e6]。不然Gurobi反馈“Unbounded objective”其实是因为目标中缺少对eta的约束绑定。正确做法是用不等式约束for k 1:N cons [cons, z_k(k) 0, z_k(k) loss_k - eta]; end cons [cons, sum(z_k)/N 1e5]; % 防止无限增大 obj ... eta (1/(1-alpha))*sum(z_k)/N;这样eta自动被推高到满足约束的最低水平不会被极端的损失样本拉爆。6.3 两阶段耦合变量在MATLAB中的传递我最初直接写了一个函数second_stage_cost(x, xi)试图在优化过程中调用。但YALMIP不支持在约束里调用任意MATLAB函数——必须把第二阶段的内层优化展开为约束。解决办法是用对偶化把第二阶段问题转化成一组线性约束加到主问题中这一步不能偷懒。如果不熟悉对偶推导建议先看几篇Wasserstein分布鲁棒的经典文献比如Mohajerin Esfahani和Kuhn的论文把对偶公式抄下来再实现。6.4 备用耦合约束容易忘记“跨时段”约束储气罐备用和机组爬坡约束是跨时段的。比如一台机组上备用预留了但下一时段如果爬坡速率不够实际无法实现额外出力。我遇到过调度结果在实时阶段不可行的情况原因就是只加了静态备用约束没加爬坡约束。解决办法是预留备用时同步加机组爬坡约束cons [cons, P_G(i,t1)-P_G(i,t) Ramp_up(i)*15 R_up(i,t)]; cons [cons, P_G(i,t)-P_G(i,t1) Ramp_dn(i)*15 R_dn(i,t)];这样才真正满足“备用可调用”的要求。6.5 计算时间优化的小技巧分布鲁棒CVaR模型需要几百个样本时求解时间可能超过十分钟。我有两个提速技巧第一压缩样本数量。不是每个时段都需要全部样本参与可以用聚类方法把相似误差模式聚成10个典型历史场景再在模糊集里用聚类中心权重计算Wasserstein球。这样精度损失不大但速度提升几倍。第二用热启动。把上一时段的决策作为MIP的初始解传入Gurobi通过ops.gurobi.Start设置初始变量值能明显减少分支次数。另外关闭对数输出verbose0能节省一点点时间但主要节省的是外部I/O。6.6 对气网动态的简化要谨慎如果项目背景是短时间内的实时调度比如5分钟间隔那天然气动态延迟不可忽略。我校验过在日内滚动调度中忽略管存效应会导致燃气机组爬坡时气网节点压力跌穿下限。所以在滚动调度版本里我额外加入了管存模型V_pipe M * dp把管道惯性等效为虚拟储气罐。不过这些扩展做起来非常绕如果只是做学术复现或毕业设计固定步长的稳态气网模型已经够用。7. 一点个人经验与扩展方向这套代码跑通之后有几个地方值得再往上走。第一可以加入碳交易机制把碳排放配额作为随机变量扩展到环境-经济联合调度。第二可以把条件风险价值换成谱风险度量比如广义扭曲风险这样能更精细地刻画尾部风险偏好。第三如果算力允许可以把电-气扩展到电-气-热综合能源系统热网惯性大备用耦合的方式会更复杂也更加考验建模功力。我个人在实际操作中最深的体会是分布鲁棒优化的价值不在于给出一个“最优解”而是在于给出一个“在多种可能世界中都不太差”的决策。这种思路非常契合电-气综合能源系统这种多主体、多时间尺度、多类型不确定源混杂的场景。如果你正打算复现相关论文或者拿它做毕设我建议先不要急着写代码把模糊集对偶推导和CVaR线性化公式手推一遍再去看代码细节。公式懂透了MATLAB实现只是时间问题。最后再分享一个小技巧调试这类模型时先把Wasserstein半径设成0先验证随机优化的调度结果是否合理再把半径一点点调大观察备用容量和成本的变换是否符合直觉。这样排错会快很多也不会被一堆耦合约束整得手足无措。希望这篇文章能帮你在电-气综合能源系统的优化道路上少踩几个坑。