做综合能源系统优化这些年我经常被同一个问题卡住风电预测误差的真实分布根本不知道凭什么直接用正态分布或者历史频率分布去做随机优化这个问题最稳妥的解法之一就是分布鲁棒优化——不要精确分布只要一批历史场景圈一个可能的分布集合在最坏分布下做决策。如果再叠加上条件风险价值去刻画尾部风险就成了标题里说的“计及条件风险价值的电-气综合能源系统能量-备用分布鲁棒优化”。这种题目在论文里往往一页公式接一页公式真正用Matlab把它从论文变成能跑的代码完全是另一件事。我最早复现这类模型时光是把“分布鲁棒CVaR电-气耦合备用优化”四件事揉进一个可求解的优化问题里就折腾了将近三周。这篇文章我把完整的建模思路、求解框架、MatlabYalmip实现要点和调试经验都写出来给打算做这个方向或者正在被论文公式折磨的同学一条能直接抄的近路。1. 问题起点风电不确定性遇上电-气耦合难点到底在哪1.1 为什么偏偏是电-气综合能源系统传统电力系统优化只盯电网这一侧但燃气机组的大量接入让电力系统的灵活性越来越依赖天然气网络。天然气管道有管容约束、有压力约束、有压缩机特性气源还有供应上限这些物理限制会直接传导到电力系统你以为燃气机组还能往上调100MW但气网在那个时段根本供不上这么多气。反过来电转气P2G设备让电力系统和天然气系统形成了双向耦合。风电大发的时候电网消纳不了P2G可以把多余电力转成天然气存进气管网或储气设施相当于给风电装了一个“缓冲池”。但这个缓冲池的容量、位置、运行策略又会影响电力系统的备用安排。所以电-气综合能源系统的能量-备用联合优化本质上是把一个“电网最优”的问题升级成“电网-气网联合最优”的问题。备用容量不再是电力系统单方面拍脑袋定的指标它必须同时受到气源出力上限、管道输送能力、储气设施状态的约束。1.2 经典随机规划和传统鲁棒优化为什么都不够用如果只用随机规划你得给出风电预测误差的精确概率分布。但实际工程里我们能拿到的往往只是有限个历史场景你猜的分布和真实分布之间一定有偏差。分布设错了优化出来的备用方案可能实际风险很高。传统鲁棒优化倒是把分布问题绕过去了但它要求所有不确定量都落在给定的集合内相当于给“最坏情况”设置了上限。这个最坏情况往往来自极端但概率极低的工况按这个标准配置备用成本会高到不切实际。分布鲁棒优化走的是中间路线假设真实分布落在以经验分布为中心的一个模糊集内然后优化“最坏可能分布”下的期望成本。它只需要历史样本又能通过模糊集半径控制保守程度所以在电-气综合能源系统这种数据有限、模型复杂的问题里特别合适。1.3 “能量-备用联合优化”到底在优化什么很多初学者把能量优化和备用优化当成两个先后执行的问题先出清能量计划再根据经验法则配置备用。这样做的弊端在于备用需求本质上来自风电不确定性而能量计划又决定了系统在什么运行点上保留备用两者必须同时决策才是全局最优。联合优化里第一阶段的决策变量同时包含机组基准出力、启停状态、上备用容量、下备用容量、气源供气计划、P2G运行计划第二阶段是在真实风电出力实现之后通过再调度去校正偏差追求最小的实时运行成本和失负荷、弃风惩罚之和。换句话说能量计划决定“系统在哪个点运行”备用容量决定“系统能在这个点附近扛多大偏差”两者放在同一个目标函数里优化才是真正的联合优化。2. 条件风险价值怎么进入优化模型从定义到线性化2.1 CVaR的定义和它在备用问题里的角色条件风险价值CVaR衡量的是“损失分布尾部”的平均水平。给定置信水平α比如α0.95VaR是5%分位点处的损失值而CVaR是损失超过VaR的那部分的期望。用公式表达就是CVaR_α(Z) min_η { η (1/(1-α)) · E[(Z - η)_] }其中 (Z - η)_ max(Z - η, 0)。为什么备用优化里特别适合用CVaR因为失负荷、弃风这类事件的特点是“平时不发生发生就是大事”。期望成本可能看起来不高但一旦风电预测误差碰上极端场景损失可能非常大。普通期望优化会轻视这种尾部风险CVaR则把它显式地放进目标函数逼着模型多留一些备用来压制尾部损失。2.2 把CVaR写进分布鲁棒目标函数的技巧我采用的模型形式是下面这种它在很多相关文献里都出现过兼顾了理论严谨性和代码可实现性min [ 第一阶段运行成本 λ·η θ ]s.t. θ ≥ sup_{P∈B_ρ(P_hat)} E_P[ Q(x,ξ) (λ/(1-α))·(Q(x,ξ) - η)_ ]其中x是一阶段决策变量ξ是风电不确定量Q(x,ξ)是第二阶段再调度成本B_ρ(P_hat)是以经验分布为中心的Wasserstein距离球λ是风险厌恶系数η是CVaR内部的辅助变量。这个形式把“期望成本”和“尾部风险成本”合并在一起θ对应最坏分布下的期望再调度成本CVaR那一项专门惩罚极端条件下的高成本。λ越大系统配置的备用越充足尾部风险越低但一阶段基准成本也会升高。2.3 辅助变量与线性化落地CVaR表达式里的( Q(x,ξ) - η )_是非线性的但引入辅助变量t就能线性化t ≥ Q(x,ξ) - η t ≥ 0于是目标里那部分变成 λ·η λ/(1-α)·t。注意t本身是ξ的函数在分布鲁棒框架下它要参与最坏分布的对偶变换不能简单当成常数。一套常用的参数我放在这里供参考置信水平α取0.95风险系数λ从0到1.0之间扫描失负荷惩罚取8000元/MWh弃风惩罚取300元/MWh。单位用标幺制更方便风功率预测误差归一化到[-1,1]区间。3. 分布鲁棒模型的核心Wasserstein模糊集与两阶段结构3.1 为什么用Wasserstein距离构造模糊集分布鲁棒优化里模糊集的构造方式不止一种比如矩模糊集、φ散度模糊集、Wasserstein距离球。我推荐Wasserstein距离球的原因是它很直观一切和历史样本“离得不远”的分布都算候选分布距离半径ρ直接控制了不确定程度。Wasserstein距离的通俗理解是要把一个概率分布P“搬运”成另一个概率分布Q平均每个质量点至少要移动多少距离。它天然考虑了场景之间的空间距离不像φ散度那样只盯着概率密度比值。模糊集写成B_ρ(P_hat) { P : W(P, P_hat) ≤ ρ }其中P_hat (1/N) Σ δ(ξ_i)是N个历史场景构成的离散经验分布ρ是模糊集半径。3.2 两阶段模型第一阶段和第二阶段的完整表达第一阶段是“决策-观望”结构里的先行者。在风电出力还没实现时要先定机组启停、基准出力、备用容量、气源出力和P2G计划。第一阶段的目标包括发电成本、备用容量成本、气源供应成本。第二阶段是在风电实际出力ξ实现后做再调度。这时可以调整机组出力来应对偏差代价是再调度成本如果调整不过来就切负荷或者弃风付出高额惩罚。第二阶段成本Q(x,ξ)就是这个再调度问题的最优值。电力侧的再调度约束包括机组出力上下限、爬坡约束、备用容量约束、功率平衡约束。天然气侧的约束包括气源出力上下限、节点气平衡、管道流量约束和储气约束。两者通过燃气轮机耗气量方程和P2G产气量方程耦合。3.3 为什么min-max-min结构必须经过重构直接看模型结构是min_x { 一阶段成本 sup_P E_P[ min_y 二阶段成本 ] }这是一个min-max-min三层嵌套问题求解器根本没法直接处理。必须先把内层的min用对偶理论换成max把问题变成min-max-max再利用Wasserstein对偶把sup_P展开成有限维形式。具体来说第二阶段如果是个线性规划它的强对偶可以写出来。假设第二阶段原始问题是min c^T y s.t. A y ≤ b B ξ, y ≥ 0对偶问题是max (b B ξ)^T λ s.t. A^T λ ≤ c, λ ≥ 0于是第二阶段的最优值Q(x,ξ)可以等价成一个关于λ的最大化问题。做完这步外层min和max-sup之间就剩下一个形式上可迭代处理的结构。这是整套代码里最关键的数学变换。4. 求解算法列与约束生成框架下的三层迭代4.1 为什么不用一步到位的单层重构严格来说分布鲁棒问题可以尝试写成一个大MILP或大LP一锤子求解但工程上往往不这么做。原因有两个第一把Wasserstein球的非线性项全部线性化会引入海量辅助变量和0-1变量问题规模爆炸第二电-气综合能源系统本身已经有大量0-1变量机组启停和分段线性约束再堆变量会让求解时间从分钟级变成小时级。实际代码里更常用的是列与约束生成CCG框架主问题只保留第一阶段决策和割平面约束子问题负责找“最坏场景”然后把最坏场景对应的约束反馈给主问题循环迭代直到上下界收敛。4.2 主问题、子问题怎么切分主问题求解的内容是min [ 一阶段成本 λ·η θ ]约束条件包括第一阶段所有可行性约束以及从子问题返回的一组割平面。每次迭代往里加一条割平面θ会越来越紧最终逼近真实的最坏分布期望成本。子问题的任务是在给定主问题解x和η的情况下对每个历史样本ξ_i求解max_ξ [ Q(x,ξ) (λ/(1-α))·t - γ·c(ξ,ξ_i) ]其中γ是Wasserstein对偶乘子c(ξ,ξ_i)是距离度量t对应CVaR尾部辅助变量。这个子问题本身还是min-max嵌套需要继续用KKT条件或者强对偶转成单层。我实际采用的做法把第二阶段LP写成对偶形式再用互补松弛条件把“原始-对偶同时可行且目标相等”写成一组带大M的线性约束最后整个子问题变成一个可以直接交给Cplex求解的MILP。对每个历史样本求解一次取最大者作为最坏场景。4.3 迭代收敛判据和容差设置CCG的上下界是这样维护的每次主问题求出一个下界因为主问题是原问题的松弛把主问题解代入子问题算出最坏场景再带回主问题重新求解前用当前x和子问题目标值计算一个上界。我用的收敛判据是|UB - LB| / |LB| ≤ 1e-3正常情况下这个模型5到15次迭代就能收敛。如果超过30次还不收敛基本可以断定是子问题里的互补松弛条件或者大M系数出了问题不要盲目加迭代次数。下面是主问题和子问题迭代的伪代码框架LB -inf; UB inf; cuts {}; for iter 1:max_iter % 主问题 [x_opt, eta_opt, theta_opt, LB] solveMaster(cuts); % 子问题对每个历史样本找最坏场景 [worse_val, worse_xi] solveSubproblem(x_opt, eta_opt); % 更新上界 UB min(UB, masterCost(x_opt) worse_val); % 添加割平面 cuts{end1} buildCut(worse_xi); % 收敛判断 if abs(UB - LB)/abs(LB) 1e-3, break; end end5. MATLABYalmip代码落地从数据准备到求解器调优5.1 算例数据准备场景生成、电-气耦合参数我用的算例是标准39节点电力测试系统加经典20节点天然气测试网络。电力系统10台常规机组其中2台替换为燃气轮机风电装机600MW用两个风电场接入天然气网络包含4个气源、2个储气设施。历史风电场景数我生成500个然后聚类缩减到50个。为什么缩减因为子问题要对每个历史场景跑一个MILP场景太多迭代一次要等太久。聚类的思路是用K-means或者场景消减算法把距离近的场景合并保留形状特征。缩减后模糊集半径需要乘一个修正系数这里的经验做法是把ρ按样本数比例做一个缩放。耦合参数是这类算例最容易出错的地方。燃气轮机的耗气量方程是二次的要分段线性化P2G的效率取60%左右即1MWh电力产出0.6MWh等价的天然气。5.2 主问题和子问题的Yalmip建模框架Yalmip建模的核心思路是把主问题当成一个带累积割平面的优化问题每轮迭代更新约束集合。下面是主问题建模的骨架% 主问题变量 u binvar(n_g, 1); % 机组启停 p0 sdpvar(n_g, 1); % 基准出力 Rup sdpvar(n_g, 1); % 上备用 Rdn sdpvar(n_g, 1); % 下备用 qs0 sdpvar(n_s, 1); % 气源基准供气 eta sdpvar(1, 1); % CVaR辅助变量 theta sdpvar(1, 1); % 分布鲁棒成本项 % 目标函数 objective sum(c_g .* p0) sum(c_rup .* Rup) sum(c_rdn .* Rdn) ... sum(c_gas .* qs0) lambda * eta theta; % 基础约束 Constraints [ sum(p0) sum(p_wind_predict) sum(p_p2g) load_electric, ... Rup 0, Rdn 0, ... p0 Rup Pmax .* u, ... p0 - Rdn Pmin .* u, ... ]; % 累积割平面 for k 1:length(cuts) Constraints [Constraints, theta cuts(k).expr]; end ops sdpsettings(solver, cplex, cplex.mip.tolerances.mipgap, 1e-4); sol optimize(Constraints, objective, ops);子问题的Yalmip建模稍微复杂一些要把第二阶段原始问题写成对偶形式并做互补松弛线性化。一个比较清晰的写法是把每个历史场景单独循环for i 1:N_scenarios xi sdpvar(n_w, 1); % 该场景的不确定变量 y sdpvar(n_y, 1); % 二阶段连续变量 lam sdpvar(n_con, 1); % 对偶乘子 z binvar(n_con, 1); % 互补松弛0-1变量 % 原始可行性 cons_sp [A * y b B * xi, y 0]; % 对偶可行性 cons_sp [cons_sp, A * lam c, lam 0]; % 互补松弛线性化 slack b B * xi - A * y; cons_sp [cons_sp, slack 0, slack M * (1 - z), lam M * z]; % 目标对偶值 CVaR辅助 Wasserstein距离惩罚 obj_sp -(b B * xi) * lam - gamma * norm(xi - xi_hist(:,i), 1); % 求解该子问题记录最坏场景值 end注意子问题目标里我取的是负号因为要把外层max转成Yalmip的min形式。这里的方向问题是初学者最容易翻车的地方建议先用一个只有单风电场、无气网的小例子验证对偶方向再扩展到完整算例。5.3 求解器参数与效率调优这类问题的计算瓶颈通常在子问题的MILP求解上。我实测下来合理设置求解器参数能让总时间下降一半以上Cplex的MIP gap设到1e-4就够了再小意义不大。开启并行计算Threads设成CPU核数减1。LP求解器选dual simplex因为子问题每次求解的热启动效果很好。另外强烈建议给子问题设置求解时间上限。比如单个子问题超过60秒就直接接受当前可行解。因为CCG迭代的容差是1e-3个别子问题求解不精确不会明显影响最终结果但会让总时间爆炸。6. 结果解读与常见调试陷阱6.1 算例结果里能挖出什么结论我用上面这套流程跑完算例后最关心的三张图是总成本随模糊集半径ρ的变化、备用容量随风险系数λ的变化、CVaR值和期望成本之间的权衡曲线。先说模糊集半径ρρ0退化成经验分布下的随机规划总成本最低但尾部风险最高ρ逐渐增大总成本上升但CVaR值明显下降说明系统用一部分经济性换取了更低的极端风险。这个曲线是审稿人和导师最爱的图也是选ρ的重要依据。风险系数λ的变化更有意思λ从0增到1上备用容量明显增加下备用的增幅小一些燃气机组备用容量的增幅受到气网供气约束的限制。做灵敏度分析时能清楚看到“电-气耦合约束成为备用瓶颈”的现象这正好说明了联合建模的价值。6.2 大M系数取值和数值稳定性互补松弛线性化里的大M是调试过程中最大的坑。M太小会切掉最优解M太大会让数值稳定性急剧恶化。我一开始图省事统一取M1e6结果求解器频繁报数值警告。做了几轮实验后总结出一个可靠做法把模型中涉及的量全部归一化到标幺值范围比如功率基准取100MW成本基准取相应数量级然后大M取100到1000之间。如果你的模型必须保留原始单位大M不能统一给一个值必须分约束给对应的量级比如气网流量的M和失负荷惩罚的M完全不在一个数量级。如果Cplex报“infeasible”或者Yalmip报“Numerical issues”优先检查大M相关的行把M缩小两个数量级再试。6.3 场景数小时模糊集半径的修正问题很多复现者会把历史场景直接从论文里抄过来但场景数量不同时模糊集半径不能直接套用。分布鲁棒优化的理论给出半径ρ和样本数N有近似关系在同样置信度下样本越少需要的模糊集半径越大因为经验分布的不确定性更高。我自己的经验是如果你的有效场景数只有50个而参考论文用了500个半径大概要放大1.5到2倍才能得到相似保守度。具体做法是跑不同ρ下的结果画成本-风险曲线结合你想承担的风险水平来定别迷信论文里的一个固定数字。6.4 燃气轮机备用受气网约束时的线性化细节燃气轮机二次耗气特性如果不做分段线性化Yalmip会把子问题变成MIQP甚至直接报错。我的做法是取3到5段线性逼近横坐标是燃气轮机出力纵坐标是耗气量分段点均匀分布在出力上下限之间。分段数太少误差大太多则增加0-1变量3到5段是性价比最高的区间。还有一个容易被忽略的约束燃气轮机提供上调备用时它预留的那部分容量对应的额外耗气也要在气网约束里预留出来。也就是说气源供气计划里要包含“备用所需气量”否则实时阶段无法兑现备用。我第一次建模时漏了这条约束结果算出来的备用方案在子问题里大量切负荷前几轮迭代上下界一直横跳就是这个原因。7. 扩展思路和实际体会这类模型搭好之后扩展方向非常多。想加储能就把储能连续变量和充放电约束并进第一阶段想把P2G的爬坡约束加进子问题也只是在第二阶段约束里多几行想从单断面扩展到24小时日前调度需要把机组启停变量变成时段维问题规模会涨不少但求解框架不用动。我个人在这个项目上最深的体会是这套代码的难点不在数学推导而在“把每一层嵌套都落实到Yalmip的变量和约束上”。对偶方向、互补松弛线性化、模糊集半径标定、割平面怎么往回填任何一个细节没对齐整个迭代就会出问题。所以我强烈建议你在完整算例之前先搭一个最小可跑通的小系统一个火电机组、一个风电场、一条天然气管道、一个燃气机组跑通CCG迭代再慢慢放大。最后分享一个小技巧把每一步的上下界保存在数组里每次迭代后画出来。UB和LB曲线如果像剪刀一样匀速收敛说明算法一切正常如果出现锯齿或者长时间不降直接去查子问题的最坏场景是不是被大M限制了。这条经验帮我省过太多无谓的调试时间。