分布鲁棒机会约束与N-1准则下的低碳经济调度复现指南

我前阵子接了个复现任务标题是“考虑N-1准则的分布鲁棒机会约束低碳经济调度”。光看名字就知道这不是个普通的经济调度安全校核、分布鲁棒、机会约束、低碳机制全搅在一起。我原本以为论文有公开代码结果只找到一篇EI收录的正文和几个模糊的公式剩下全靠自己补。断断续续折腾了两周总算把Matlab代码的主流程跑通也让结果过了蒙特卡洛检验。这篇博文不是论文翻译而是我从零复现时踩过的坑、想通的逻辑、以及最后怎么确认“这个结果真的对了”。如果你也是在复现类似DRCC调度模型或者想理解N-1准则怎么和分布鲁棒机会约束共存可以参考一下。1. 先别急着写代码这个题目到底在算什么1.1 低碳经济调度把碳排放进“成本账”传统经济调度追求发电成本最低一般只包含煤耗成本、购电成本、启停成本、弃风弃光惩罚这几项。低碳经济调度则多做了一件事把碳排放也变成成本的一部分。具体做法并不唯一常见的是碳税还有配额-交易机制以及近些年论文里很爱用的阶梯碳交易机制。不管哪种形式最后结果都会改变最优出力分布——高碳机组会被压低出力低碳机组甚至储能会被多调用风电光伏的消纳优先级也会提升。复现时最容易忽略的点是“碳排放不只在目标函数里出现”。有的模型会给系统一个碳排放上限这就变成一个全局约束有的模型把碳排放配额分配到机组变成每台机组的约束还有的模型引入碳捕集装置让机组发电的同时产生一部分可交易的碳配额。它们在代码里的写法完全不同。我建议你先把论文的目标函数和约束逐条列出来搞清楚碳到底出现在哪几个位置再动笔写代码。1.2 N-1准则的作用不只是“多一条约束”N-1准则并不新鲜电力系统规划和安全校核里一直在用。简单说就是系统里任意一个元件退出运行后电网仍然不能出现负荷损失、线路过载、电压越限这些问题。在确定性模型里做N-1校核就是把每个预想事故单独跑一遍潮流这很容易理解。可一旦放到随机调度里事情就变复杂了。新能源出力、负荷预测误差本来就是不确定的事故状态又是一个离散的随机事件。如果一个约束既要覆盖连续随机量又要覆盖各种N-1事故那它就不再是普通的确定性潮流约束了。我复现时最大的感受是N-1在确定性模型里是“枚举重算”在DRCC模型里是“把离散事故状态和连续随机变量一起塞进机会约束”。这比单纯加一条备用容量约束要难一个数量级。1.3 分布鲁棒机会约束不确定性建模的折中为了说清楚分布鲁棒机会约束我习惯把几种做法放一起看。建模方式对不确定性的假设典型后果随机规划已知精确概率分布分布不准时结果可能严重失稳传统鲁棒优化只知道不确定量的边界最坏情形下覆盖成本偏高分布鲁棒优化真实分布落在一个“模糊集”内未知但受限保证最坏分布下满足约束比鲁棒便宜机会约束允许小概率违约风险预算可控但分布不准确时评价会失真DRCC就是分布鲁棒优化和机会约束的结合我先圈定一个模糊集把所有“看起来合理”的分布都放进去然后要求在这个模糊集内的最坏分布下安全约束仍然能以不低于1-ε的概率成立。说白了你不需要精确知道风电预测误差的真实概率但你要为一大类可能的概率分布做好准备。1.4 复现这个工作需要的知识包如果你只想抄一段代码跑通那确实不用懂全部推导。但你想修改边界条件、换测试算例、解释结果至少得准备这几块凸优化和对偶理论、混合整数线性规划、概率论基础、Matlab的矩阵化编程习惯。尤其是Wasserstein距离和对偶转化它们是DRCC的核心跳过这层直接看代码很容易在模糊集半径或者机会约束的方向上栽跟头。2. 分布鲁棒机会约束先把模糊集和机会约束的转化搞明白2.1 为什么直接样本分布靠不住很多第一次接触DRCC的人会问我已经有历史场景了直接估计风电误差分布然后写一个普通的机会约束不就行了吗理论上可以但实际样本量通常不够。你要是只有几百个历史场景估计出的分布尾部误差很大而安全约束恰恰关心尾部。风光预测误差的厚尾特性又比较明显一旦把尾部低估算出来的调度方案在真实场景里违反安全约束的概率会明显超过设定值。分布鲁棒思想的出发点就在这里不把样本分布当作真实分布而是假设真实分布位于“以样本分布为中心的一个邻域”里。这个邻域就是模糊集。只要这个模糊集半径给得合理无论真实分布怎么在邻域内波动机会约束的违约率都能得到控制。这就是为什么DRCC比纯随机规划更能让人放心。2.2 Wasserstein距离球是首选模糊集模糊集的构造方式很多有矩模糊集、KL散度球、Wasserstein距离球。我做这个复现时最常用的是Wasserstein距离球因为它在场景数据上很自然而且对偶之后能转化成可以交给优化求解器处理的形式。Wasserstein距离可以通俗理解为“搬土成本”它衡量两个分布之间最小的运输代价。以样本经验分布为中心半径ρ内的所有分布构成模糊集。当历史场景数N增加时ρ会按大约O(N^{-1/2})的速度缩小。这就解释了为什么DRCC最终会向随机规划收敛数据越多模糊集越小越敢依赖样本分布数据越少模糊集越大越偏向保守。2.3 机会约束等价转化的思路机会约束一般写成这种形式Pr{ g(x, ξ) ≤ 0 } ≥ 1 - ε其中x是调度决策ξ是随机变量。如果直接把ξ的所有可能分布都考虑进来这个问题无法直接求解。但把真实分布限定到Wasserstein模糊集之后原约束可以通过对偶转化为一组有限维约束。转化后会出现辅助变量和0-1变量模型变成混合整数线性规划或混合整数二次规划。我复现时没有去手推每一个对偶步骤但有一件事必须搞清楚转化后的约束中那个ε到底出现在不等号的哪一边以及ρ是加在约束的哪一项上。论文里符号稍微不一致代码就会差出一个数量级。这个只能靠对照原文逐行确认没有捷径。2.4 参数到底怎么选两个参数最关键一个是允许的违约概率ε另一个是模糊集半径ρ。ε一般取0.05或0.10太小会把模型搞得过于保守太大又失去了安全意义。ρ的取值需要根据场景数量和置信水平来定常见的经验公式带有sqrt(d * log(1/β) / N)这种结构其中d是随机变量维度β是置信度参数。实际操作时我更建议做一个灵敏度扫描固定ε让ρ从0开始逐步增大观察目标函数是否单调变差。如果目标值完全不变可能是约束没有激活如果变动太剧烈可能说明随机变量单位不一致。我当时就是吃了单位的亏——风电误差用的是MW约束方程里用的是标幺值导致ρ无论取多大都像没取一样。3. N-1安全校核与机会约束是怎么合流的3.1 N-1预想事故集怎么来做N-1校核第一步是确定预想事故集。如果枚举所有线路和机组组合数量会非常大而且很多元件退出运行对系统影响很小。合理做法是先跑一遍确定性潮流筛掉那些即使断线也不会让任何线路越限的元件保留可能影响输电断面、可能引起负荷损失的N-1事故。在代码实现里事故集通常是一组“元件索引”。比如线路故障就记录这条线路的两端节点和电抗机组故障就记录这台机组所在节点和容量。不要小看这一步事故集不准确后面的机会约束再完美也没用。3.2 故障状态下的潮流方程怎么写正常状态下的直流潮流可以写成P_line B * (θ_i - θ_j)线路断开后这条线路不再承担潮流同时系统拓扑发生了改变。简单粗暴的方式是重新构造节点-支路关联矩阵再重新算B矩阵。但调度优化里不能让矩阵在每个场景里动态变化所以更多是用大M法把正常状态和故障状态的约束写在一起。我习惯引入一个0-1变量u_f表示事故f是否发生。正常时u_f0约束启用原有线路方程故障时u_f1把线路潮流强制为0同时用修改后的节点功率平衡方程代替原来的平衡方程。写成代码时本质上是在原始约束右侧加上M * u_f这一项让故障时原约束自动松弛掉。这个大M的取值很关键取太小会切掉可行域取太大则数值不稳定。3.3 把事故状态放进机会约束这一步是整个复现里最难的地方。连续随机量是风电和负荷预测误差离散随机量就是哪些元件发生N-1故障。两者同时出现在机会约束里时数学上会变成一个混合随机规划。我看过的几类论文里处理方式有两种倾向一种是把N-1事故单独拿出来在确定性层面枚举只把连续随机量放进机会约束另一种是把离散事故状态也作为随机变量让机会约束同时覆盖故障事件。前者实现起来更容易后者更贴合“统一不确定性”的表达但需要额外处理离散分布与模糊集的耦合。你复现的时候一定要先搞清楚参考论文用的是哪一种因为这两种写法的最优成本和代码结构差别非常大。我自己先按第二种思路建了一个模型结果因为事故状态和连续随机变量的联合分布不好处理求解器频繁报错后来改成“离散枚举连续DRCC”才把主流程跑通。3.4 别把每个N-1事故都单独立一个机会约束一个容易犯的错误是为每一条预想事故分别写一个Pr ≥ 1-ε的机会约束。这样模型倒是清晰但保守性会被严重放大。因为所有事故都要求同一个很高的保证概率相当于要求系统在最恶劣事故组合下也几乎不越限成本会向传统鲁棒优化靠拢。更合理的做法通常是给所有预想事故一个共享的风险预算或者让不同事故的违约概率根据发生的严重程度分别设定。实现时可能需要引入“联合机会约束”的近似处理比如Bonferroni近似把总风险预算拆分到各个事故约束上。这样写出来的模型既保留了N-1安全要求又不会让成本爆炸。4. 低碳目标怎么建模碳交易、碳捕集和成本线性化4.1 碳交易机制的三种做法低碳经济调度里碳成本建模大概有三种。第一种是固定碳税每吨碳排放直接乘一个单价第二种是配额-交易机制系统先获得免费碳排放配额实际排放超过配额就要买低于配额可以卖第三种是阶梯碳交易超额越多单价越高形成分段递增的碳价曲线。从数学角度固定碳税最简单是个线性项配额交易需要加入买卖变量阶梯碳交易则需要分段线性化会引入整数变量和辅助连续变量。论文里用阶梯碳交易的不少因为它的减排激励更强而且能体现“越超额越贵”的政策逻辑。复现时不要急着把所有细节堆上去我建议先用固定碳税把模型跑通再切换成配额或阶梯机制这样出错了也容易定位。4.2 碳捕集与“净排放”如果模型里包含碳捕集装置事情会更复杂一步。碳捕集装置在捕集二氧化碳的同时要消耗电能所以机组发出的电分为两部分一部分供给系统负荷一部分用于捕集。发电毛出力减去捕集能耗才是送入电网的净出力。同样碳排放也分总排放和净排放总排放由发电毛出力决定捕集掉的碳可以从总排放里扣减扣完之后才是参与碳交易的净排放。这个双重“净”概念特别容易在代码里写错。我见过好几个初版实现要么在功率平衡里忘了扣捕集能耗要么在碳成本里直接用总排放参与交易最后结果比实际高出一大截。建议你把每个机组单独建一个数据结构字段里同时写上P_gross、P_ccs、P_net、E_total、E_captured、E_net在代码里显式传递这些量不要靠公式隐式换算。4.3 机组成本二次曲线分段线性化机组耗量特性通常是一条二次曲线燃料成本是出力的凸函数。DRCC模型本身已经带了0-1变量和大量辅助变量如果目标函数再直接保留二次项求解压力会很大。实际复现时我一般把二次曲线分段线性化把机组出力范围切成若干段每段用一个线性斜率近似。分段线性化的代码实现并不复杂。每台机组引入分段出力变量P_seg和一个或几个0-1变量目标函数累加“段出力乘以该段斜率”同时保证各段按顺序填充。分段数太少近似误差大分段数太多求解时间成倍增加。我通常取3到5段足够说明方法趋势又不至于把算例卡死。4.4 目标函数的量纲和权重把电能量成本、弃风弃光惩罚、碳交易成本加在一起时最容易被忽略的是量纲。电价通常用元/MW·h碳排放价格通常用元/吨。机组排放强度单位是吨/MW·h所以要先把排放量折算成与电价相同的能量成本单位再进入目标函数。我建议你先把不考虑碳约束的确定性经济调度跑一遍看各项成本大概在什么量级。比如总发电成本是100万元碳成本只有几百元那后面无论怎么调碳价影响都很小如果碳成本高到几百万说明碳价设置不合理低碳目标已经完全压过了供电成本。调参时先固定一台机组的典型排放强度算一笔直观账再整体微调碳价系数。5. Matlab代码复现的关键步骤从搭数据到调通5.1 数据准备阶段不要用Excel硬填复现这种模型第一步是准备算例数据。节点系统一般包含节点编号、支路两端编号、电抗、容量、机组成本系数、负荷水平、风电出力序列等。不要把所有数据都塞进Excel再手工读取最好在你自己的脚本里用矩阵和结构体统一管理。我在代码里喜欢把所有基础数据放进几个结构体bus、branch、gen、scenario。scenario里存的是历史风电场景和误差场景。构建约束时尽量用矩阵运算替代循环。比如节点功率平衡可以写成A_bus * P - L 0而不是一行行累加。矩阵运算不仅快而且不容易出现索引错位。5.2 模型搭建顺序我习惯按这个顺序搭模型先声明关键变量再写目标函数然后逐类添加约束。约束的顺序一般是功率平衡、机组出力上下限、爬坡约束、线路潮流约束、N-1事故约束、碳约束、机会约束。每加一类约束后就检查一次变量个数和矩阵维度不要等到所有约束写完了再一次性调试。对了机会约束这一步要确认哪些约束需要被“机会化”。如果连普通功率平衡都被写成了机会约束模型会变得无比难解而且物理上也不太合理。一般情况下被机会化的是线路潮流安全约束和备用容量约束功率平衡仍然保持确定性或场景化约束。下面是当时主循环的伪代码结构方便参考% 声明关键变量 P_g sdpvar(n_gen, 1); % 常规机组出力 P_w sdpvar(n_wind, 1); % 风电计划出力 s_s sdpvar(n_slack, 1); % 辅助变量 u_f binvar(n_con, 1); % N-1事故状态变量 % 目标函数 Objective fuel_cost(P_g) carbon_cost(E_net) curtail_penalty; % 基础约束 Constraints []; Constraints [Constraints; sum(P_g) sum(P_w) - L 0]; Constraints [Constraints; P_g_min P_g P_g_max]; % 线路潮流约束正常态 Constraints [Constraints; -Cap Bf * theta Cap]; % N-1事故约束 for k 1:n_con Constraints [Constraints; branch_fault_constraint(k, u_f, M)]; end % 机会约束转化后的约束 Constraints [Constraints; drcc_safety_constraint(epsilon, rho, scenarios)];代码里的drcc_safety_constraint就是一个函数它把Wasserstein模糊集下的机会约束等价形式写进去。建议先让这个函数返回空约束跑通其他部分再逐步加入真正的DRCC约束。5.3 模型求解的工程化处理模型转化后是MILP节点规模稍大时求解时间会明显上升。工程上有几个办法缓解一是给求解器设置时间上限和相对间隙比如允许1%的次优性二是先求解一个“松弛LP”把0-1变量固定成连续变量看看可行性三是按事故集分批求解先看有没有明显不可行。如果求解器报模型不可行不要急着怀疑理论。先做两件事检查所有大M约束是否在正常状态下松弛再检查辅助变量的上下界是否正确。很多时候不可行只是因为某个辅助变量没有给定上下界或某条支路容量写成了负数。5.4 我在复现时候踩过的具体坑第一大M取值。N-1事故约束里故障时要让原线路潮流约束自动失效需要在约束右侧加大M。各支路潮流范围差很多不能用一个全局大M。我后来按每条支路容量上限的若干倍取M并做了多组敏感性测试。第二随机变量单位。风电出力误差是MW约束里如果用了标幺值模糊集半径ρ也要换算成标幺值。否则ρ变成纯数字看起来调了很多实际效果约等于零。第三ε的方向。有的论文写Pr{违约}≤ε有的写Pr{满足}≥1-ε移植公式时多一个负号结果完全相反。我自己就在一个约束里把1-ε写成了ε1导致模型极度保守最后检查了半天才发现。第四故障线路的关联矩阵。线路断开后功率平衡必须用修改后的关联矩阵而不是原矩阵。如果不改相当于故障线路仍然连着节点N-1约束形同虚设。我当时漏了这一步任何故障都不会导致越限结果自然“好得离谱”。5.5 验证代码是否正确的三步法模型写完之后目标值再好看也不能直接相信。强烈建议做三步验证。第一步是极限验证。把机会约束里的ε设到接近1比如0.99也就是几乎允许所有场景违约模型应当退化成普通确定性模型。如果结果和确定性模型差很多说明某个约束写重了或者没正确松弛。第二步是蒙特卡洛检验。生成几千个随机场景把求出的最优调度代进去统计安全约束的真实违约率。如果真实违约率明显高于设定的ε说明模糊集半径太小或转化写错了如果远低于ε说明模型偏保守。第三步是参数灵敏度。固定其他条件把模糊集半径ρ从零逐渐加大看总成本是否单调上升。分布鲁棒模型天然有这个性质如果没有基本可以断言约束没有被正确激活。6. 结果应该长什么样用“逻辑链”判断复现对错6.1 成本构成表怎么看复现完成后第一件事不是看总成本而是看成本构成是否合理。我以一个标准30节点系统为例列出典型趋势注意具体数值因算例而不同不要拿这里的数直接当基准。模型总成本常规机组成本碳交易成本弃风率确定性模型低低高低随机规划较低较低中等中分布鲁棒机会约束中中低中传统鲁棒优化高高低高这个表格反映的是逻辑趋势越保守的模型总成本越高常规机组被调用得更多弃风率更高但因为低碳机组出力更多碳交易成本反而会降低。如果你复现出的DRCC结果比随机规划还便宜那一定有问题。6.2 结果排序的底层逻辑分布鲁棒机会约束的定位很明确它应该比随机规划保守比传统鲁棒优化乐观。如果最优成本排不出这个顺序大概率是下面几种原因模糊集半径太小导致DRCC退化成普通随机规划机会约束的违约概率设得过大约束基本不起作用或者模糊集半径太大DRCC直接撸成了鲁棒优化。我习惯把结果按“确定、随机、DRCC、鲁棒”四个模型各跑一遍画成柱状图。看到DRCC的柱子确实夹在随机和鲁棒中间心里才踏实。6.3 检查解是否真的满足N-1目标值和蒙特卡洛检验都过了还要专门验证N-1。做法很简单取DRCC求出的最优机组出力和风电出力人为断开预想故障元件重新跑一次潮流计算检查所有支路负载率和节点电压是否越限。如果任意一个预想故障下出现过载说明N-1事故约束没有正确生效。这一步尤其重要。因为机会约束允许小概率违约N-1安全可能只是“大概率安全”。你需要在结果报告里写明在哪些场景组合下允许违约违约概率是多少。这样才算真正复现了一个有工程意义的调度模型而不是只交出一个目标值。7. 复现这类论文我最后想说的两件事7.1 把论文里的符号统一写一页复现之前先花半天把论文里的所有符号整理到一页纸上。这个工作非常枯燥但特别有效。DRCC模型里符号特别多ε、ρ、N、K、ξ、x、u等等论文里同一个字母在不同章节可能指不同含义。一旦摘抄代码时理解错一个符号后面所有结果都会偏。我通常是先建一个符号表把“变量的含义、维度、单位、在代码中的变量名”四列对应起来。代码里尽量沿用这个命名不要一会儿用P_g一会儿用Pg否则调试到后半夜会谢自己。7.2 留好每一步的过程数据这种复现工作几小时跑一次算例很常见。我建议把每一步的中间变量、目标值、约束违约率都保存到文件里而不是只打印在命令行。这样切换参数时可以快速对比也方便回溯是哪一步改坏了结果。如果你也打算复现这个题目我可以分享的经验是不要期待一次跑通把模型拆成模块先确定性、后随机、再分布鲁棒每个阶段都配一个验证脚本。DRCC真正难的不是公式而是把公式变成你不会怀疑的代码。等蒙特卡洛检验通过、排序逻辑对上的那一刻你会觉得前面踩的坑都值了。