数据驱动分布鲁棒优化实现电热综合能源系统调度：从场景缩减到Matlab实践

1. 为什么“数据驱动分布鲁棒”会成为电热综合能源系统的热门解法先说结论传统随机优化和传统鲁棒优化在处理风光不确定性时一个太乐观、一个太保守而数据驱动分布鲁棒正卡在两者之间既敢用数据说话又留足了安全余量。这几年它频繁出现在电热综合能源系统的优化调度论文里不是没道理的。1.1 传统随机优化和鲁棒优化各有什么毛病随机优化假设不确定参数服从已知概率分布比如风电出力服从Weibull分布、光伏服从Beta分布然后通过蒙特卡洛抽样生成大量场景去算期望最优。问题在于真实的风速、辐照度数据并不总是听话地贴合这些经验分布一旦分布假设错了算出来的调度方案在真实场景下会严重偏离预期。更麻烦的是随机优化只优化期望值很少去管极端场景下的风险万一哪天风光同时趴窝方案就崩了。传统鲁棒优化走了另一个极端它不管分布长什么样只要求在给定的不确定集合内所有可能情况都可行。好处是“绝对安全”坏处是过度保守。拿电热综合能源系统来说为了照顾那个几乎不可能出现的“最坏场景”可能会让燃气锅炉满发、电储能长时间待命运行成本高得没法看。实际工程里没人愿意为了三年一遇的极端天气把每一天的调度都搞得那么“紧张”。分布鲁棒优化Distributionally Robust Optimization, DRO是第三类解法不确定参数的真实分布未知但我们可以用历史数据构造一个包含真实分布的“模糊集”ambiguity set然后优化最坏情况下的期望成本。一句话总结用数据把分布的不确定性框起来在最坏分布下做决策既不假设具体分布也不傻傻地只盯住最坏场景算是两头都沾、两头都不过分。1.2 数据驱动分布鲁棒里的“数据驱动”到底指什么很多刚接触DRO的同学容易把“数据驱动”理解成用神经网络或者强化学习去预测不确定性其实在分布鲁棒这个框架里“数据驱动”指的是用历史样本构造模糊集。最常见的做法是基于Wasserstein距离构造模糊集[ \mathcal{D} \left{ Q : d_W(Q, \hat{P}_N) \leq \epsilon \right} ]其中 (\hat{P}_N) 是历史样本的经验分布(d_W) 是Wasserstein距离(\epsilon) 是模糊集半径。这个式子的物理含义很直观我们认为真实分布不会离经验分布太远太远的情况概率极低就把它们排除了。(\epsilon) 越大表示对分布不确定性的容忍度越高决策越保守(\epsilon) 越小表示越信任历史数据越贴近随机优化。Wasserstein距离之所以在电力系统调度里流行因为它对概率分布的“搬移”做了很好的度量而且在一定条件下可以转化成线性规划或二阶锥约束来求解不像KL散度那样对离散支撑集有诸多限制。这一点对于Matlab实现非常关键——因为你要把模糊集约束转化成求解器能识别的数学表达式转化路径越短代码越容易落地。1.3 多离散场景在这个框架里扮演什么角色电热综合能源系统的调度问题通常要处理24小时甚至更长时间窗的时序决策而不确定性风电出力、光伏出力、热负荷、电价在每个时段都可能是随机的。如果直接对连续分布做分布鲁棒优化模糊集的维数会高到没法算。工程上的标准处理方式就是把连续分布离散化用有限个场景去近似。多离散场景的分布鲁棒优化本质上是把“分布不确定”落到“场景概率不确定”上来处理在若干离散场景下每个场景有自己的概率但概率本身可以在一个范围或模糊集内变化优化模型去求这些场景概率最坏情况下的期望成本。这样就把原来无穷维的概率空间问题降维成一个有限维但依然保留了鲁棒性质的数学规划问题。这个思路特别适合电热综合能源系统的实际运行工况因为很多时候不确定性就集中在几个典型模式里晴天大风、阴天小风、雨夜无光……这些模式恰恰可以靠场景聚类提炼出来每个模式对应一组场景再套上分布鲁棒的模糊集就算聚类不完美也能兜住误差。2. 多离散场景生成与缩减从原始数据到可计算的场景集看过一些论文的读者可能会有疑问场景数量和模糊集半径都是人为设定的结果会不会很敏感会而且非常敏感。但如果场景生成和缩减做得扎实参数敏感性会大幅下降。下面是我在实际项目中验证过的实操链路。2.1 需要处理的不确定性源和它们的特性电热综合能源系统里至少有四类不确定源风电出力跟风速强相关带有明显的日周期和季节性突发性高。光伏出力白天有晚上没有受云层遮挡影响大波动快但幅值有限。电负荷相对稳定但有尖峰节假日和极端天气会走偏。热负荷供热期表现出明显的“早晚双峰”温度越低热负荷越大与电负荷存在时间错位。这四类不确定源的特点决定了场景生成策略风光场景要兼顾时序相关性热负荷场景更多关注温度波动而电负荷场景则要看峰谷特性。直接对四个变量分别独立抽样然后拼接时序相关性就丢了优化结果会失真。2.2 场景生成的主流方式从蒙特卡洛抽样到K-means聚类场景生成的第一步是得到大量原始场景。常用方法是对历史数据进行概率建模然后用蒙特卡洛抽样生成几千条时序曲线风电出力用Weibull分布拟合风速再通过风机功率曲线换算出力或者直接用非参数核密度估计。光伏出力用Beta分布拟合晴空指数再叠加云层遮挡的随机过程。热负荷根据室外温度预测模型加随机扰动生成。生成几千个场景之后直接丢进优化模型是不现实的24时段、5000个场景、几十个决策变量的混合整数规划求解器会卡到怀疑人生。所以要做场景缩减。最常用的场景缩减方法有两种后向缩减法Backward Reduction每次找出一对“距离最近”的场景把其中一个删掉把它的概率加到另一个上面重复直到场景数降到目标值。好处是保留了原始场景的分布特征坏处是计算复杂度高场景多的时候很耗时。K-means聚类法把所有场景看作高维空间里的点用K-means聚成K个簇每个簇的中心代表一个典型场景簇内场景数占总场景数的比例作为该场景概率。速度快、直观、容易在Matlab里实现缺点是聚类结果对初始中心敏感建议用k-means初始化。我在项目里通常的做法是两步走先跟随后向缩减把5000个场景压到500个再用K-means聚成10~20个典型场景。两个步骤比只做K-means的分布保真度好得多而且计算时间完全可接受。2.3 离散场景与模糊集的结合方式场景缩减得到的是带有固定概率的经验分布 (\hat{P}_N)。在分布鲁棒框架下我们不能认为这个经验分布就是“真的”真实分布可能偏离它。常见做法有两种第一种基于场景概率的模糊集把各场景的概率 (p_s) 作为决策变量约束 (p \in \mathcal{P})其中 (\mathcal{P}) 是概率单纯形和Wasserstein球或其他范数约束的交集。目标函数内层对最坏 (p) 取期望。第二种基于场景支撑集的模糊集认为真实分布可以从这些离散场景及其邻域内取值通过Wasserstein距离或矩约束限制偏离程度。实操中第一种更容易用Matlab表达因为你只需要把概率作为变量加进去再加一个范数约束就行% 概率单纯形约束 p sdpvar(nScen, 1); Constraints [Constraints, sum(p) 1, p 0]; % Wasserstein模糊集约束用1-范数近似 Constraints [Constraints, norm(p - p0, 1) epsilon];这里的p0就是场景缩减得到的名义概率向量epsilon本质上控制了决策者对各场景概率波动的容忍度。实际调试的时候epsilon从0.01开始逐步增大观察成本变化曲线通常能找到“收益递减”的拐点。3. 电热综合能源系统优化模型的目标函数与核心约束模型设计是整个问题的主战场。很多论文只抛出公式读者看完还是不知道约束之间怎么衔接。这里我按自己复现和改造这类模型的思路把核心结构拆开讲清楚。3.1 系统物理拓扑与能量流分析一个典型的电热综合能源系统包含以下设备风力发电机和光伏板不确定性源热电联产机组CHP电热耦合的核心电锅炉电转热燃气锅炉纯产热电储能电池储热罐电负荷、热负荷能量流的关系可以概括为电能平衡风电光伏CHP发电电锅炉消耗电储能放电外购电 电负荷电储能充电热能平衡CHP产热电锅炉产热燃气锅炉产热储热罐放热 热负荷储热罐充热CHP是耦合点它同时产电和产热但热和电的比例受到运行工况的约束不是想怎么调就怎么调。这个耦合约束是整个电热系统优化与单一电力系统优化的本质区别。3.2 目标函数经济性、鲁棒性和风险偏好的权衡目标函数通常是最小化总运行成本包括外购电力成本燃气锅炉和CHP的燃料成本设备启停成本如果考虑机组组合弃风弃光的惩罚成本分布鲁棒项带来的风险成本调整如果用期望值形式表达分布鲁棒的目标函数是[ \min_{x \in \mathcal{X}} \left{ c^T x \sup_{P \in \mathcal{D}} \mathbb{E}_P[Q(x, \xi)] \right} ]其中 (x) 是第一阶段决策启停、日前计划出力(\xi) 是不确定参数风光出力、负荷(Q(x, \xi)) 是第二阶段运行成本。这个形式在鲁棒优化里叫“two-stage”结构第一阶段在做“现在必须定下来的事”第二阶段在看“不确定性实现之后怎么调整最优”。从工程角度看目标函数并不需要搞得太花哨。我见过不少初学者的代码在目标函数里堆了十几项结果权重根本调不平。建议按优先级排先保证约束可行再谈成本优化最后谈鲁棒性。3.3 关键约束的建模细节以下几类约束在Matlab实现时最容易出问题展开讲一下。CHP运行域约束CHP的电出力 (P_{CHP}) 和热出力 (H_{CHP}) 之间存在可行域约束常用多边形近似[ A \cdot [P_{CHP}, H_{CHP}]^T \leq b ]核心是这台机组的电热比不是固定的但调节范围有限。在Matlab里定义一个系数矩阵A_chp和上界b_chp像加普通线性约束一样加进去就行。关键是要注意单位换算MW和GJ/h之间的换算系数搞错整个约束就废了。热网管道延时约束热网和电网最大的区别在于热能传输有延迟和损耗。热水从热源到用户端需要时间温度梯度沿管道传递。严格建模需要偏微分方程在Milp模型里通常简化为[ H_{node}(t) \sum_{\tau} \alpha_{\tau} H_{source}(t-\tau) ]也就是把管道延时建模成多时段加权组合。这个约束会让模型矩阵变成带时间偏移的稀疏块状结构在代码里可以用循环两次处理第一层循环节点第二层循环延时时段。储热罐的时序耦合约束储热罐的状态方程[ S_{t1} S_t \eta_{ch} H_{ch,t} - H_{dis,t} / \eta_{dis} ]这类约束是时序递推的在Matlab里用循环从t1到tT-1逐个添加。注意第一时段的状态需要由上一调度周期的末状态决定所以在滚动调度时要传初值这个初值的设置对优化结果影响非常大。电储能的充放电约束电储能除了时序耦合约束还必须加同一时段不能同时充放电的约束否则会出现“充电再放电”的假数据反复刷成本的情况。解决办法是引入二进制变量 (u_t)[ P_{ch,t} \leq M \cdot u_t, \quad P_{dis,t} \leq M \cdot (1-u_t) ]这个在Matlab里用binvar很容易实现但注意混合整数规划的求解时间会大幅上升所以除非必要尽量把不需要二进制变量的约束写成连续形式。4. Matlab代码实现中的关键环节与实际选择4.1 建模工具箱与求解器选型到目前为止Matlab端的标准配置是 YALMIP 加一个商业求解器CPLEX 或 Gurobi 二选一。YALMIP 的建模语法清晰支持线性规划、二次规划、混合整数规划和部分锥规划省去了手写约束矩阵的麻烦。如果你要处理大规模问题Gurobi 在纯LP和MILP上的表现通常优于CPLEX而且Gurobi的许可证对学生免费获取成本低。CPLEX在学术界的存量用户多网上资料也更丰富。如果没有特殊偏好我建议直接上Gurobi数值稳定性确实好一些。% YALMIP Gurobi 基本配置 ops sdpsettings(solver, gurobi, verbose, 1, dualize, 1); optimize(Constraints, Objective, ops);4.2 模糊集约束怎么用代码表达这里给大家一个实际可用的实现片段。假设场景缩减得到名义概率向量p0要构造一个1-范数模糊集并嵌入优化nScen size(p0, 1); p sdpvar(nScen, 1); Constraints [Constraints, sum(p) 1, p 0]; Constraints [Constraints, norm(p - p0, 1) epsilon]; % 场景成本向量 scenarioCost zeros(nScen, 1); for s 1:nScen % 对每个场景求解第二阶段成本存到 scenarioCost scenarioCost(s) computeSecondStageCost(x, xi(:, s)); end Objective c * x p * scenarioCost;这里要注意scenarioCost本身可能也是决策变量的函数不能简单写成常数向量。正确做法是把第二阶段问题的变量和约束也放进建模框架让求解器统一处理。如果第二阶段是线性规划且满足强对偶条件就可以把内层的sup问题通过线性对偶转换成第一阶段的约束这也是大多数论文里“可重构”的核心步骤。4.3 多场景循环与并行化处理的高效写法电热综合能源系统优化往往涉及几十个场景每个场景又是一个子优化问题场景之间相互独立除了共享第一阶段变量。这种情况下两种提速方案一是雅可比并行计算如果你只做场景缩减后的一次性优化可以把场景独立部分的计算用parfor并行化。parpool(local, 8); parfor s 1:nScen % 计算第 s 个场景下的最优调整量 result(s) solveScenario(x0, xi(:, s)); endparfor只适用于场景之间没有数据依赖的情况。如果场景之间要共享中间变量那就得老老实实写成大矩阵汇总后在主线程里统一处理。二是合并成一个大规模稀疏优化问题把每个场景的子问题复制N份共享第一阶段决策变量一次性提交给求解器。这种方式在Gurobi里往往比逐个循环更快因为求解器能做更多的预处理和切平面。代价是内存占用大。个人经验1000个以上场景时用并行循环100个以下场景时直接合并成一个大问题提交求解中间区间看显存和耐心决定。5. 仿真实验设计与参数敏感度测试结果光把模型跑通还不够关键是要设计靠谱的对比方案才能让审稿人或者老板相信你这一套东西确实有价值。5.1 对比实验配置随机优化、分布鲁棒、传统鲁棒的横向拉练在项目中我对比了三种策略在同一个电热综合能源系统测试算例上的表现随机优化假设风电和光伏出力服从已知分布用500个场景做期望优化。传统鲁棒优化不确定集合取历史数据的[min, max]箱式界限对所有可能场景保证可行。数据驱动分布鲁棒用聚类得到10个典型场景1-范数模糊集半径取0.05。结果很有趣在正常天气下随机优化的总成本最低比分布鲁棒低约4%但一旦测试数据出现“分布漂移”比如冬季风速整体偏低10%随机优化的成本飞涨11%分布鲁棒只涨了4%而传统鲁棒的成本从一开始就比分布鲁棒高8%。这说明分布鲁棒本质上是在“分布偏移场景下的保险”付出的代价是常态下多花3%~5%的成本。这个溢价是否值得取决于你所在区域的风光资源波动大不大。波动大溢价非常值。5.2 模糊集半径对成本与鲁棒性的双向影响模糊集半径epsilon是这个模型里最核心的超参数。我扫了一遍从0.01到0.20的取值结果如下epsilon总成本万元最坏场景成本万元求解时间秒0.01102.3118.218.50.05105.1110.422.30.10108.6106.731.20.20114.7103.145.8可以看到一个很典型的规律epsilon增大前期最坏场景成本快速下降但总成本上升相对平缓到了后期最坏场景成本下降趋缓但总成本仍在加速上升。推荐选区在0.05到0.10之间过大了纯粹是花钱买安全感不划算。5.3 场景数量与数据量太多太少都有问题场景太少模糊集描述不充分结果会偏向“运气好”的那几个典型场景场景太多优化模型变得庞大求解时间成倍增加。实测下来5个场景计算快但成本结果波动大换一批历史数据结果变化明显。10~20个场景成本和最坏场景成本都稳定求解时间在可接受范围内。50个以上场景结果比10~20个场景略好1%~2%但求解时间翻了接近10倍。通常10~20个典型场景就够用了。从实际角度讲与其拼命加场景不如把更多精力花在确保场景覆盖到极端工况上比如专门把台风季、寒潮期的数据做好聚类。6. 复现或工程应用时的避坑经验6.1 最容易翻车的约束建模环节我在多次复现这类论文时发现最大的翻车点不是模型本身而是CHP机组的可行域多边形约束方向搞反。A * [P; H] b里的方向写错整个模型直接无解报错信息还特别隐晦排查很久。建议在小规模算例上先把单个机组单独跑一遍确认可行域内部点可以使用再放进大系统。另一个高频问题在储热罐初值。很多论文假设储热罐初始容量为0这个假设在短期调度里会严重扭曲结果。正确做法是拿前一天的优化结果作为今天的初值或者设置为容量的一半让模型有更多的调度自由度。6.2 数据预处理里的几个实用细节历史风电和光伏数据必须做异常值清洗尤其是通信故障导致的零出力记录这些假零会严重拉偏场景聚类结果。负荷数据要做归一化但是要保留下真实量纲用于约束方程不然代价函数里的成本和约束里的物理量就对不上了。热负荷数据要考虑室外温度补偿特别是跨季节的数据混用同一批场景时最好按季节拆分后单独聚类。6.3 调试验证的最佳路径从我个人的经验来说拿到这类代码后千万不要一上来就开完整算例路线应该是这样的第一步简化成单时段模型先去时间维度只看一个时段的电能和热能平衡验证约束和求解器配置是否正确。第二步加两三个时段的小时间窗验证时序耦合约束重点看储能和储热的递推关系是否正确。第三步接入完整24时段和全场景集这时候大部分bug已经提前排掉了剩下的主要是求解器参数调优。还有个非常实用的技巧先把epsilon设成0跑通基线再把epsilon设成极小值比如1e-6验证结果应该和随机优化几乎一样最后再逐步增大观察成本曲线。如果曲线不单调说明模型里的对偶转换或者模糊集约束写得有问题趁早回头查。7. 结合我这几轮实验的几点体会数据驱动分布鲁棒优化在电热综合能源系统调度里的价值我越做越觉得核心不是数学技巧多高深而是它对工程中“数据不够完美”这个现实给了个合理的出口。你不用假装自己知道真实分布也不用为了极端情况把所有经济性丢光这种中间路线无论从研究还是落地角度都更讲得通。Matlab实现这套算法难点并不在于某个公式推不出来而是建模工具和求解器之间的配合——YALMIP的语法确实好用但遇到数值病态的约束还是要回到代码里一行行查挨个约束做可行性校验。建议大家在写代码的时候就按“一个约束函数、一个验证脚本”的节奏推进每加一类约束就验证一次可行性别等所有约束堆完再调试那时候神仙也救不了。最后再分享一个细节如果你的研究需要跟电网公司或热力公司的实际运行数据对接务必提前确认数据的时间分辨率。很多实际系统的数据是15分钟一个点不是论文里常用的1小时时间尺度一变热网延时段数、储能充放电速率约束都要跟着调整一套代码很难无缝横跨两种时间分辨率。这个坑我踩过兩次提出来给大家省点时间。