分布式鲁棒优化与联合机会约束的电力调度MATLAB实现

我得先给这个标题祛个魅。分布式鲁棒优化、联合机会约束、能量与储备联合调度这三个词叠在一起乍一看像是又一个高不可攀的电力系统优化论文题目但本质上它解决的是一个非常实际的问题当前天的风电预测曲线出来后第二天机组怎么开、出力定多少、留多少备用容量才能在风光出力波动和预测偏差面前既保证系统安全、又不让成本高到离谱。我最近在实验室的MATLAB测试环境里完整跑了一遍这个流程把分布式鲁棒优化DRO的模糊集建模、联合机会约束的转化、以及能量和备用的联合决策揉进了一个可复现的框架。这篇文章就把我的建模思路、MATLAB实现过程和踩过的坑一次性说清楚。如果你正在做新能源消纳、调频备用优化或者想找一个机会约束规划的落地案例这篇内容应该能帮你省掉几周试错时间。1. 为什么能量与储备必须放在同一个优化里1.1 联合调度的核心价值传统的调度模式是先做能量调度也就是经济调度让机组出力满足负荷预测然后再根据预测偏差人工估算备用需求通常取最大单机容量或者负荷的某个百分比。这种两步走的做法在新能源占比低的时候问题不大但一旦风电光伏大量接入预测误差和波动性大幅上升人工定备用的方式就开始失灵。我举一个实际算例里的场景某时段预测风电出力300MW如果只做了能量调度机组按这个值安排出力但实际风电出力可能只有250MW这时50MW的缺额需要实时平衡市场高价购买一次两次还能接受整个调度周期内频繁出现就会显著推高运行成本。反过来如果备用定得过多全天都在为用不上的容量支付成本同样是浪费。能量与储备联合调度的核心价值就是在一个优化问题里同时决定机组出力、上备用容量和下备用容量让它们共享机组容量、爬坡速率这些物理约束做到在能量计划里预设备用空间在备用决策里反哺能量分配。这样备用水平不再拍脑袋定而是根据不确定性的大小和成本结构自动权衡。1.2 三种不确定性处理路线的取舍研究这个课题时最先要回答的问题是怎么处理风电预测误差这个不确定性。我对比过三条技术路线这里直接说结论。第一种是随机规划SP假设预测误差服从某个已知分布比如正态分布然后用场景法采样求解。这条路线的问题在于真实预测误差的分布往往比正态分布更厚尾而且均值、方差都在随时间变化所谓已知分布其实建立在错误假设上。一旦真实分布和假设分布偏差较大求出来的解就可能低估风险。第二种是鲁棒优化RO不假设分布只要误差落在一个不确定集合内模型就保证所有集合内的情况都安全。这个方法看起来保险但实际计算时集合边缘的极端场景概率极低几乎不会发生强行照顾这些场景的结果就是备用容量被抬高成本明显偏大。我在同一个算例里测过纯RO模型的备用成本比DRO模型平均高出大约15%到20%。第三种就是分布式鲁棒优化DRO。它的思路是我不假设精确分布但我假设真实分布落在某个模糊集里。这个模糊集可以由历史数据的均值、协方差、支撑集等统计信息来刻画。优化目标变成在模糊集内最坏情况分布下的期望成本最小化。打个生活类比SP像是一个自称很懂本地天气的专家你完全按他报的概率决策RO像是不管天气预报直接按全年最冷的那天准备衣服DRO则像是你知道本地气候的大致范围按这个范围内最坏的可能性准备但不考虑外星天气。三者的信息需求、保守度和计算复杂度跨度很明显实际选型时我最终选了DRO原因是风电预测误差的分布很难精确建模但历史数据足够用来估计统计量DRO刚好能用上这部分信息又不至于过度保守。1.3 为什么在MATLAB环境里做我选择在MATLAB测试环境里实现这套策略不是因为它比Python高级而是因为它确实适合这类优化模型的快速迭代。YALMIP这个建模工具箱加上Gurobi或者Mosek这样的商业求解器几行代码就能描述一个带二阶锥约束或者半定约束的优化问题这在Python里需要手动组织稀疏矩阵工作量完全不是一个量级。另一个原因是调试体验。求解器报错之后我直接在命令行里检查决策变量、约束残差、对偶变量几秒钟就能定位问题。在快速验证模糊集半径取多少合适这类参数敏感性问题时这种交互式工作流很顺手。另外学术界公开的电力系统算例数据很多都是.m格式的像Matpower工具包里的case30、case118拿过来就能直接配合使用。MATLAB环境还有一个容易被忽略的好处代码可以无缝对接Simulink仿真模型。我后续要验证调度策略在动态过程中的表现就可以把优化结果直接喂给Simulink里的发电机、负荷和风电模型形成完整的闭环测试链路。2. 联合机会约束的本质与转化2.1 联合约束比单约束强在哪里很多初学者第一次接触机会约束时会混淆两个概念。单个机会约束处理的是某一个不等式约束成立的概率要大于某个阈值比如功率平衡约束成立的概率要大于0.95。但实际调度中系统安全不是靠单个约束保证的而是多个约束同时满足才能成立。举个例子假设系统里有10个安全约束如果每个约束单独满足的概率都是0.95在它们相互独立的理想情况下所有约束同时满足的概率是0.95的10次方约等于0.598远远低于预期的0.95。真实电力系统里约束之间是耦合的违规还会连带影响其他约束实际联合违规概率比独立假设下可能更高。联合机会约束就是把这些约束打包要求它们同时成立的概率不低于一个阈值。这才是真正意义上的系统级安全。我在建模时把功率平衡约束、备用容量约束和线路潮流约束放在同一个联合机会约束里确保它们只在很低的概率下被违反而不是各自为政。2.2 从概率约束到可求解的确定性形式联合机会约束直接求解几乎不可能因为它涉及高维概率积分所以我把它转换成可以交给求解器处理的确定性或凸约束。主流方法有三种我分别说一下取舍。第一种是Bonferroni近似把联合概率约束拆成多个单机会约束把总违规概率分配到各个约束上然后对每个单约束单独处理。这个方法实现简单但结果保守只适合约束数量很少的情况。我在早期版本里试过因为约束多每条约束分到的概率很小导致备用容量被抬得很高。第二种是采样近似也叫场景法。用大量随机场景把机会约束转成在每个场景下都成立的确定性约束然后通过场景削减技术减少计算量。这个方法在随机规划里很常见但在DRO框架下需要仔细处理模糊集和采样之间的关系场景太多计算量爆炸场景太少概率保证又不够。我在算例阶段使用它来做验证不作为主求解方式。第三种是CVaR近似。CVaR叫条件风险价值含义是最坏的那一小部分情况下的平均损失。用它来近似机会约束的好处是机会约束要求的是一段概率区间都满足而CVaR把偶尔违规一次变成违规的平均程度要小这比布尔式的概率判断更平滑数学性质也更友好。在DRO框架下CVaR近似可以和应用模糊集的期望操作自然地结合起来最终转化成一个凸优化问题。我最终采用的就是这个方案。转化后的约束形式大致是这样的把所有关于不确定变量的约束写成g(x,ξ)≤0的形式然后把联合机会约束替换成一个包含正部函数[g(x,ξ)]_的CVaR约束再加上一个辅助变量来线性化正部函数。这个过程中会引入一些松弛变量和线性化约束但整体模型仍然保持凸性适合用内点法求解。2.3 模糊集的构建与参数调节模糊集的构建是DRO模型的核心也是最容易出问题的部分。我使用的是矩模糊集它利用历史数据估计预测误差的均值向量和协方差矩阵然后用这两个统计量来约束模糊集内所有可能的分布。更严谨地说模糊集要求任意分布与经验矩的偏差在一定范围内这个范围用参数gamma来控制。gamma取0时模型退化成随机优化认为矩信息完全准确gamma越大模糊集越大越保守求解结果越倾向于配置更多备用。我在项目里花了很大功夫调这个参数。用过一组典型数据测试gamma从0.2逐步增大到2.0时总成本上升曲线先陡后缓在gamma等于1.0到1.5之间时输出最稳定越过1.5后成本上升明显但安全性改善不大这就是过度保守的信号了。除了矩模糊集还可以用Wasserstein距离模糊集它是以经验分布为球心、半径为半径的球。这种模糊集在处理高维问题时表现不错但需要调距离半径而且转化成确定性模型时会引入更多辅助变量。我在项目中做了一组对比发现矩模糊集在计算效率和保守度之间更平衡最后还是用了它。需要注意支撑集也不能忽略。即使模糊集定义了矩信息理论上仍可能存在概率质量落在物理不可能区域的分布所以要在模糊集里显式加入支撑集约束把预测误差限制在合理范围内。3. MATLAB环境下的完整实现流程3.1 工具箱配置与求解器选型实际跑代码之前先把环境配好。我用的是MATLAB R2022b配合YALMIP做建模求解器用的是Gurobi 10.0。安装完后在MATLAB里执行yalmiptest检查一下确认Gurobi能被正确识别。我的调用配置是这样写的ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.TimeLimit, 600, gurobi.MIPGap, 0.0001);选Gurobi是因为它处理二次目标函数、二阶锥约束和混合整数问题的能力都很稳而且在YALMIP里的接口很成熟。如果你的环境里没有GurobiMosek也可以SDPT3则更适合小规模问题但不推荐在大规模DRO模型里用它。另外记得把YALMIP安装到非系统盘的纯英文路径下我之前放在D盘一个含中文名的文件夹里导致工具箱加载报错排查了半天才发现是路径问题。3.2 算例系统与数据准备验证环境用了修改版的IEEE 30节点系统基础数据来自Matpower的case30.m然后做了修改在节点14和节点24接上两个风电场额定容量分别是200MW和150MW总负荷水平取原始系统的1.2倍模拟高比例新能源场景。机组参数我整理成了一张表写代码时直接加载。这里列出了6台机组的关键参数。机组容量(MW)a($/MWh)b($/MWh)爬坡率(MW/h)备用成本($/MW)G12000.003520503.0G21000.004022403.2G3800.004525303.5G4500.005028204.0G5300.006030154.5G6200.007035105.0风电数据这块我用一个ARMA模型生成500个预测误差场景代表历史样本然后用k-means聚类削减到30个代表性场景。这30个场景用来估计均值向量和协方差矩阵也就是模糊集的输入。负荷曲线取典型日的24点数据由基本负荷加峰荷系数构成。所有数据统一用MW为基本单位成本用美元避免单位和量纲混乱。3.3 核心建模代码框架下面这段代码是我实际运行的核心框架我简化掉了部分约束细节保留主干结构。它体现了DRO模型的三个层次决策变量定义、确定性约束、机会约束转化后的约束。% 系统参数定义 nGen 6; T 24; nWind 2; % 决策变量机组出力P、上备用Ru、下备用Rd P sdpvar(nGen, T, full); Ru sdpvar(nGen, T, full); Rd sdpvar(nGen, T, full); % 风险变量和辅助变量 theta sdpvar(1, 1); s sdpvar(nConstraint, nScen, full); % 目标函数发电成本备用成本DRO最坏期望调整成本 objective sum(sum(a.*P.^2 b.*P c)) ... sum(sum(ruCost.*Ru rdCost.*Rd)) ... theta; Constraints []; % 确定性约束机组出力上下限、爬坡约束、备用上下限 Constraints [Constraints, Pmin P Pmax]; Constraints [Constraints, -ramp P(:,2:T)-P(:,1:T-1) ramp]; Constraints [Constraints, 0 Ru RuMax, 0 Rd RdMax]; % 功率平衡基准场景约束 % baseLoad是典型日负荷baseWind是预测风电 Constraints [Constraints, sum(P Ru - Rd, 1) baseLoad - baseWind]; % 联合机会约束的CVaR近似转化 % mu和Sigma来自历史场景估计gamma是模糊集半径 % 这里用对偶变量引入模糊集的矩约束 Constraints [Constraints, mu*linVar gamma*sqrt(Sigma*quadVar) ... theta]; % 求解配置 ops sdpsettings(solver, gurobi, verbose, 2, gurobi.TimeLimit, 600); optimize(Constraints, objective, ops);注意到我上面把__mu__和__Sigma__的表达式做了简化。实际建模中YALMIP里表示SIGMA的二次项需要用到矩阵变量和二次约束比如用aux变量和约束__aux (Sigma * z) * inv(Sigma) * (Sigma * z)__来线性化。表述方式虽然细节多但核心思想是把DRO对最坏分布的期望成本upper bound转化为几个对偶变量和统计量的确定性表达式这些表达式本身是凸的。3.4 参数调优与结果验证模型跑通后我花了两天时间做参数扫描。置信水平alpha从0.01、0.05、0.1这组里测试发现alpha取0.05时成本和安全性最平衡。gamma按0.5、1.0、1.5、2.0扫描验证了前文说的过了1.5以后收益递减。验证机会约束是否真的满足我用的是蒙特卡洛法生成1000个独立测试场景把优化得到的调度方案代入统计违规率。实测结果是设定5%置信水平时实际违规率大约在3%到4%略微偏保守这是CVaR近似的正常表现。如果违规率低于设定值太多说明模型过于保守该减小gamma或alpha。我把不同置信水平下模型的表现整理成了一张表方便后续对比。置信水平总成本($)实际违规率求解时间(s)0.01128400.8%1280.05125103.5%1100.10122608.2%960.201198017.5%82这组数据说明置信水平越高成本越高违规率越低符合预期。但0.20那组实际违规率17.5%距离设定值很近工程上已经不安全了不建议采用。4. 常见问题与排查技巧实录4.1 求解器报Infeasible的排查流程第一次跑完整模型时Gurobi直接返回infeasible整整浪费了我两天时间。后来我总结了一套排查流程现在遇到这个问题都是按这个顺序走。先关掉机会约束相关的所有约束只保留确定性约束看模型能不能求解。如果连确定性模型都不可行问题出在基础约束里重点检查功率平衡方向和备用容量方向是不是搞反了。我犯过的错误是把上备用Ru的方向写反导致功率平衡约束在基准场景下就无解。如果确定性模型可行再加入机会约束用YALMIP的__analyze__命令查看约束类型确认没有出现非线性项。还有一个常见原因是模糊集参数gamma设得过大导致模糊集里包含了无法满足支撑集的分布可行域为空。把gamma调小到0.5再试通常能定位问题。4.2 数值病态与单位不统一电力系统优化模型容易得数值病态根源经常是单位不统一。我一开始混用了MW和W目标函数里既有0.001量级的成本系数又有1e6量级的变量求解器直接罢工。后来统一成MW和美元并给目标函数除以一个基准值10000做了归一化发现求解稳定性大幅提升。另一个经验是约束的标度要从建模伊始就控制。比如线路潮流约束里如果出现10的6次方系数通过修改变量为单位化处理尽量把系数控制在0.1到100之间。这样求解器的预处理阶段不会丢掉有效约束也不会出现数值假不可行。4.3 求解时间爆炸的应对DRO模型的求解时间是个大问题尤其是当场景数和约束数增加的时候。我在500个原始场景上跑过一次求解时间超过30分钟完全不实用。后来通过k-means场景削减把场景数降到30个求解时间压缩到2分钟以内同时通过对协方差矩阵的对角化简化二次项处理进一步加速。如果业务上需要实时性或频繁求解可以考虑用外逼近算法迭代求解。先解一个松弛版本然后根据最优解不断增加割平面。这个思路在学术文献里叫L-shaped方法用在DRO模型上同样有效只是实现复杂度高一些。我在项目里用到的是YALMIP的默认求解流程加上场景削减已经能满足测试需求。4.4 与确定性模型对比时的常见陷阱对比DRO模型和确定性模型的优劣不是简单地看成本数值。如果只用确定性模型的那一套成本来比较会得出确定性模型更经济的错误结论因为它没有考虑违约带来的风险成本。我的做法是同时统计三个指标总成本、实际违约率、以及违约引起的调整费用。只有三个指标放在一起看才能反映模型的实际价值。另外要保证对比时用的是同一套数据包括同一组风电场景、同一组负荷曲线和相同的机组参数。我在早期对比时重用了一次不同风功率的场景集导致结果偏差特别大后来统一数据源才得到合理结论。4.5 排错速查表我把常见坑整理成了一个速查表方便你对照排查。现象可能原因解决办法求解器报Infeasible备用方向反了或模糊集过大关闭约束逐步排查调小gamma求解时间超长场景数过多、变量规模过大场景削减对协方差做对角化数值病态单位不统一统一MW目标函数归一化结果过于保守gamma偏大或复杂性低减小gamma增大alphaCVaR近似导致违约率偏低近似模型天然偏保守校准后再降低置信水平最后分享两点个人经验这套DRO加联合机会约束的框架跑通后我最大的体会是不要一上来就追求复杂的模糊集。先把最简单的矩模糊集和Bonferroni分解跑通理解每一部分的作用再逐步换成更精准的Wasserstein模糊集和CVaR近似。我一开始直接上了完整的CVaR近似加联合机会约束出了数值问题完全不知道是哪一环的锅后来拆开重新搭建才真正搞明白。另外多花时间在数据预处理上。很多收敛问题和数值病态根源不在模型而在数据没有做单位统一、没有做归一化、没有做场景削减。这些工作枯燥但省下的是后面排查的一个个通宵。最后分享一个调试小技巧在YALMIP里求解后用__check__命令查看每条约束的残差能快速定位是哪条约束在作妖。配合__assign__函数手动赋一组可行解可以确认模型的逻辑有没有写错。我靠这个方法把一段一直报Infeasible的模型在十分钟内定位到了备用容量约束的大小写颠倒这个经验应该也适用于你。