最近几年在电力系统优化方向只要你投过稿或者审过稿一定会频繁撞见分布鲁棒优化这个词。如果再叠加数据驱动电热综合能源系统这几个标签基本就是当下最热的一类题目了。我自己的体会是这个方向之所以火不是因为它数学上有多炫而是它确确实实解决了传统随机优化和鲁棒优化在工程应用里的两难问题既不过分依赖概率分布假设又不会保守到决策没法用。这篇博文就从一个能直接跑起来的Matlab代码方案入手把数据驱动多离散场景分布鲁棒的建模逻辑、两阶段求解框架、CG迭代实现以及算例设计讲透适合正在做电热综合能源系统优化、或者想把DRO方法落地到代码里的研究生和工程师参考。我默认你具备基础的优化建模知识懂一点YALMIP知道什么是两阶段随机规划。如果这些还不熟建议先跑通一个简单的两阶段随机规划再回来看本文不然节奏会有点快。1. 电热综合能源系统优化的核心痛点为什么传统随机优化不够用很多人一开始接触电热综合能源系统优化时第一反应是列约束、写目标函数、调求解器但忽略了最麻烦的一个问题风电出力和热负荷的不确定性到底怎么刻画。这个问题的答案直接决定了优化结果能不能在真实环境里落地。1.1 源荷双侧不确定性带来的维数灾难电热综合能源系统的不确定性不只来自风电、光伏这些电源侧还来自热负荷和电负荷的用户侧行为。热负荷尤其麻烦它的波动跟天气、建筑热惯性、用户作息都有关很难用一条简单的概率曲线描述。如果要做精细的日内调度每个时段都要考虑不确定性一个24时段的问题每个时段取10个离散场景组合起来就是10的24次方量级的场景树直接求解是不现实的。所以工程上通常做一个简化用场景削减技术把成千上万的历史样本聚成少量典型场景每个场景代表一种可能的不确定性实现。这种做法在随机规划里很成熟但问题是——当你把场景聚出来后还是需要给每个场景配一个概率而这个概率本身就是从历史数据里估计出来的估计误差是不可避免的。1.2 随机优化与鲁棒优化的两难困境传统的随机规划Stochastic Programming假设场景概率是精确已知的然后最小化期望成本。这个假设在理论上是清晰的但实际中你拿到的历史数据永远是有限样本估计出来的概率分布跟真实分布之间总有偏差。更糟的是当样本数不够时优化结果会对概率估计很敏感出现在历史数据上表现很好、在真实运行中却不尽如人意的过拟合现象。鲁棒优化走的是另一个极端它要求所有不确定性在给定集合内的最坏情况下都可行。这个思路很稳健但代价是决策过于保守。对电热综合能源系统来说热负荷的波动范围如果取一个很大的区间优化出来的结果会让系统长期处于高成本、低经济性的状态这在工程上很难接受。提示分布鲁棒优化的本质是在这两者之间找一个中间点——允许概率分布在某个模糊集内变动但要求决策在这个模糊集内的最坏概率分布下也有良好表现。1.3 分布鲁棒优化的数学直觉一个类比我经常用一个比喻来解释分布鲁棒优化随机优化是你坚信天气预报说下雨概率70%就按70%带伞鲁棒优化是不管天气预报说什么按百分之百要下暴雨来准备分布鲁棒优化则是天气预报说70%但你承认这个数字可能有误差于是假设真实概率在40%到90%之间都算数然后按最坏的那种情况来带伞。这个思路放在电热综合能源系统里具体操作就是用历史数据构造一个包含多个可能概率分布的集合模糊集然后优化目标函数在这个模糊集上的最坏期望成本。因为模糊集是数据驱动的样本越多、质量越高模糊集就可以取得越小决策也越精准——这就是数据驱动四个字的含义。2. 数据驱动多离散场景模糊集的构造逻辑模糊集Ambiguity Set是整个分布鲁棒优化模型的心脏。它的构造方式直接决定了模型的保守程度、求解难度和实际效果。这篇博文聚焦的多离散场景路线是工程上最实用的一类构造方式。2.1 从历史数据到离散场景K-means聚类与场景削减第一步是把历史数据变成可计算的离散场景。假设你有过去一年的风电出力和热负荷历史数据按小时采样就是8760组样本。直接用这么多样本做优化不现实所以要做场景削减。我常用的做法是先做K-means聚类把8760个样本聚成K个典型场景每个场景的概率初始设为该簇样本数占总样本数的比例。这里K的取值需要权衡K太小场景代表性不足模糊集的基准分布偏差大K太大求解速度会明显下降。从我的测试经验看对于24时段、单节点电热系统K取20到50之间比较合适如果系统规模大可以适当减小。以下是K-means聚类生成离散场景的Matlab代码片段我用的是自带函数方便复现% data: T x N 矩阵T为时段数N为历史样本天数 % 这里把每天的时序曲线当作一个样本做时序聚类 rng(42); K 30; % 场景数 [idx, C] kmeans(data, K, MaxIter, 500, Replicates, 10); % C: T x K每个列向量是一个典型场景曲线 p0 histcounts(idx, (1:K1)-0.5) / length(idx); % 经验分布概率注意如果数据里同时含风电、热负荷、电负荷三类不确定性聚类时要放在同一个特征空间里做联合聚类而不是分别聚完再拼凑否则会破坏场景之间的相关性结构。这一点我在早期踩过坑分开聚类得到的场景集在代入优化模型后约束满足率明显偏低。2.2 模糊集的三要素支撑集合、基准分布与距离度量数据驱动的分布鲁棒模糊集一般由三个要素构成第一是支撑集合Ξ也就是所有可能场景的集合。在多离散场景框架下支撑集合就是聚类得到的K个典型场景。第二是基准分布p₀通常取历史样本频率分布。它是模糊集的中心。第三是距离度量用来定义哪些分布算在模糊集内。最常用的是1-范数或∞-范数约束形式如下D₁ { p ≥ 0, Σp_k 1, ‖p - p₀‖₁ ≤ θ₁ } D∞ { p ≥ 0, Σp_k 1, ‖p - p₀‖∞ ≤ θ∞ }其中θ₁和θ∞是模糊集半径度量了对基准分布的允许偏离程度。这个半径不是随便拍的它应该随样本量N增大而减小理论上最优收敛速率大致是O(1/√N)量级。实际调参时可以从一个初始值开始按比例扫描看不同半径下结果的变化趋势。2.3 矩信息模糊集与Wasserstein模糊集的取舍除了离散概率范数模糊集学术界还常用两类矩信息模糊集和Wasserstein距离模糊集。矩信息模糊集约束分布的均值和协方差在一个估计区间内数学形式漂亮但对实际数据来说一阶矩和二阶矩往往不足以刻画风电出力的多峰分布特性容易丢失分布形状信息。Wasserstein距离模糊集是近几年最火的方向它以经验分布为球心、Wasserstein距离为半径画一个分布球。好处是可以直接从数据构造且有一些理论上的有限样本保证。但缺点是Wasserstein球的约束往往把问题变成更复杂的半无限规划虽然可以用对偶转化但Matlab实现难度和工作量都更高。相比之下多离散场景概率范数模糊集的结构最简单可以直接写成线性约束与两阶段混合整数规划模型嵌合时不会改变问题类别用CG算法求解时子问题的对偶形式也非常干净。这也是我推荐用它作为入门方向的原因。提示如果你的审稿人或者导师特别看重理论深度可以考虑在离散场景基础上加矩约束比如Σp_k ξ_k 的期望落在给定区间内这样既保留了线性结构又提升了理论说服力。2.4 多离散场景的离散到底指什么这里值得单独解释一下。多离散场景分布鲁棒这个说法在不同论文里有不同含义但主流理解是不确定性参数的支撑集合是有限离散点集而概率分布在这K个离散点上是未知的、属于某个模糊集的。也就是说随机性有两个层次——外层是哪个场景会发生内层是某个场景发生时不确定参数取什么值。当每个离散场景本身还是一个连续区间时问题就变成了离散连续混合支撑的DRO那求解复杂度会显著上升。本文讨论的方式是把每个聚类中心直接当作该场景的代表值即场景内部不再有波动。这在工程上是一种常见近似。如果你需要更精细可以对每个聚类簇内部再做一层box不确定集但那样子问题的结构就从LP变成鲁棒LP迭代求解时计算量增加不少。我建议先把基础版本跑通再考虑扩展。3. 电热综合能源系统的优化模型建立有了模糊集接下来要把电热综合能源系统的运行优化问题写成两阶段分布鲁棒模型。这里我以一个小型园区级系统为例含一台热电联产机组CHP、一台燃气锅炉、一台电锅炉、一个储热罐外加从上级电网购电。这个配置麻雀虽小五脏俱全能覆盖电热耦合的核心特性。3.1 设备建模CHP、电锅炉、储热罐的运行约束CHP机组是电热系统的核心耦合设备。它的电出力和热出力之间存在可行域约束我用一个简化的线性四边形可行域来表达避免引入非线性项P_chp_min ≤ P_chp(t) ≤ P_chp_max H_chp_min ≤ H_chp(t) ≤ H_chp_maxH_chp(t) ≤ α₁·P_chp(t) β₁ H_chp(t) ≥ α₂·P_chp(t) β₂这些约束刻画了CHP热电比可变但受限的物理特性。电锅炉EB负责把电能转化成热能模型相对简单H_eb(t) η_eb · P_eb(t) 0 ≤ P_eb(t) ≤ P_eb_max储热罐用一阶能量平衡方程描述S_HS(t1) S_HS(t) η_ch · H_charge(t) - H_discharge(t)/η_dis - L_HS(t) 0 ≤ S_HS(t) ≤ S_HS_max其中L_HS(t)是储热损失η_ch和η_dis分别是充放热效率。注意储热罐的充放热不能同时进行这需要引入二进制变量这一项会让模型变成MILP。3.2 目标函数与电热网络平衡约束目标函数是所有设备在一个调度周期内的总运行成本。第一阶段的成本包括CHP燃料成本、购电成本第二阶段再调度成本包括调整出力产生的额外费用、弃风惩罚和热负荷削减惩罚。两阶段分布鲁棒模型的紧凑形式如下minₓ cᵀx max_{p∈D} Σ_k p_k · Q(x, ξ_k) s.t. Ax ≤ b其中x代表日前决策变量机组启停、出力计划、储热罐充放热计划Q(x, ξ_k)是场景k下的第二阶段最优再调度成本其定义为Q(x, ξ_k) min_y dᵀy s.t. Wy ≥ h_k - T_k x, y ≥ 0这里的ξ_k代表聚类得到的第k个场景包含风电出力和热负荷的数值。等式平衡约束包括电功率平衡和热功率平衡P_chp(t) P_eb(t) P_wind(t, ξ) P_buy(t) P_load(t) H_chp(t) H_eb(t) H_boiler(t) H_discharge(t) - H_charge(t) H_load(t, ξ)注意第一阶段的决策变量会同时出现在两个平衡方程中这正是电热耦合的体现。3.3 两阶段模型的决策变量划分一个问题哪些变量放第一阶段哪些放第二阶段我的经验是机组启停、热电比模式选择、储热罐的充放热计划这类需要提前一天确定的变量放在第一阶段而风电实际出力与预测值偏差引起的出力调整、弃风量、热负荷削减量放在第二阶段。这样划分符合电力系统日前调度实时调整的运行机制。需要特别提醒的是储热罐的充放热状态二进制变量放第一阶段会显著增加主问题的整数变量数量。如果求解速度不理想可以考虑把充放热状态松弛为一个连续变量加上线性化约束代价是模型精确性稍有降低。工程上这个取舍通常是可以接受的。4. Matlab求解架构CG主问题-子问题迭代与代码实现两阶段分布鲁棒模型最常用的求解方法是列与约束生成算法CG也叫CCG核心思想是把原问题拆成主问题和子问题反复迭代每次把子问题识别出的最坏场景作为新的约束加入主问题。相比Benders分解CG在处理离散变量时收敛性要好很多迭代次数也更少。4.1 算法整体流程CG的流程可以概括为以下几步初始化选定初始场景集通常用基准分布里的所有K个场景令下界LB-∞上界UB∞迭代次数l1。求解主问题得到最优解(x^l, η^l)更新下界LB max(LB, cᵀx^l η^l)。固定x^l对每个离散场景求解第二阶段问题得到Q(x^l, ξ_k)。求解子问题即外层关于概率分布p的最大化问题得到最坏分布p*和对应的最坏期望成本F(x^l)。更新上界UB min(UB, cᵀx^l F(x^l))。如果UB - LB ≤ ε则停止并输出最优解否则把识别出的最坏分布及其对应的场景约束加入主问题l l1回到第2步。这里的关键在于第3和第4步。我的做法是先分别求出每个场景下的Q(x^l, ξ_k)这是个普通的LP可以用YALMIP逐个求解也可以批量向量化然后代入外层问题F(x^l) max_p Σ_k p_k · Q(x^l, ξ_k) s.t. Σ_k p_k 1 ‖p - p₀‖₁ ≤ θ₁ p ≥ 0这是一个只有K个变量和少量约束的线性规划求解极其快。同时根据对偶理论最坏分布p*一定落在模糊集的某个顶点上这保证了算法收敛也简化了切割平面的构造。4.2 主问题的YALMIP建模与实现主问题是带有置信约束的MILP。YALMIP实现时需要注意每个迭代轮次加入的新约束要对应一个最坏场景索引而不是把全部场景的约束一开始就全部写入否则主问题规模会随着迭代逐步膨胀。以下是主问题的核心建模片段% x: 第一阶段变量, eta: 辅助变量(epigraph) x sdpvar(n_x, 1); eta sdpvar(1, 1); % 基础约束 Ax b Constraints [A * x b]; % 每个迭代加入两个约束: % eta sum_k p^*_k * d * y_k (期望成本约束) % y_k 是与最坏场景相绑定的再调度变量 for i 1:length(worst_scenarios) k worst_scenarios(i); y sdpvar(n_y, 1); % 新再调度变量 Constraints [Constraints, ... W * y h(:, k) - T * x, ... y 0, ... d * y eta]; % 或 eta sum_k p_k * d * y 的线性化 end Objective c * x eta; Options sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); Diagnostics optimize(Constraints, Objective, Options);注意CG每次迭代会引入一组新的第二阶段变量y而不是复用旧的。这是因为每个场景的再调度决策是独立的新加入的约束需要新的变量来承载。这会让主问题规模线性增长但实际迭代次数通常很少10到20轮左右所以总体可控。4.3 子问题的max-min转化与对偶处理子问题的核心难度在于max-min结构。固定x^l后内层Q(x^l, ξ_k)本身是LP外层是对p的LP。求解顺序上先内后外是可行的因为内层K个LP相互独立可以并行求解。如果希望从子问题中提取对偶信息来构造更强的切割平面可以把内层LP写成对偶形式Q(x^l, ξ_k) max_λ λᵀ(h_k - T_k x^l) s.t. Wᵀλ ≤ d, λ ≥ 0把对偶形式代入外层问题后得到F(x^l) max_{p, λ_k} Σ_k p_k · λ_kᵀ(h_k - T_k x^l) s.t. Σ_k p_k 1, Wᵀλ_k ≤ d, λ_k ≥ 0, p ∈ D这个问题的目标函数里出现了p_k与λ_k的乘积项是双线性的。不过由于p和λ_k之间没有耦合约束且每个λ_k的对偶可行域独立可以分两步求解先用内层LP解出每个场景的最优对偶变量λ_k*再代入外层求最坏分布p*。这其实就是上面流程第3、4步的数学依据。在做敏感性分析时最坏分布p*对应的那些λ_k*就是切割平面的关键信息。标准CG切割是直接把最坏分布对应场景的约束加入主问题这是最简单可靠的实现方式。4.4 数据驱动部分的关键代码段数据驱动部分把模糊集约束写成YALMIP可识别的形式。以1-范数模糊集为例% p: 概率分布变量 (K维), p0: 经验分布 p sdpvar(K, 1); theta1 0.05; % 模糊集半径需要调参 Constraints_p [sum(p) 1, p 0, norm(p - p0, 1) theta1]; % 外层最坏分布求解: F max_p sum(p .* Qk) F -inf; Qk_vec zeros(K, 1); % 每个场景下的第二段成本 for k 1:K Qk_vec(k) solve_second_stage(x_l, xi(:, k)); % 内层LP end ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints_p, -Qk_vec * p, ops); % 目标取负转为最小化 F value(Qk_vec * p); p_star value(p);这段代码执行起来非常快K50时一般不到0.1秒。整个CG迭代的总时间主要花在主问题MILP求解上所以如果你想提速重点应该放在减少主问题整数变量规模上。5. 算例设计怎么构建让人信服的对比实验代码跑通只是第一步论文和项目报告里真正有说服力的是算例设计和对比实验。这一节我分享一套我验证过、可以直接套用的实验方案。5.1 测试系统与数据准备我建议用改进的IEEE 33节点配电网6节点热网耦合系统作为标准测试平台。电网上挂CHP机组、电锅炉和风电场热网上有燃气锅炉和储热罐通过CHP和电锅炉实现电热耦合。这个系统规模适中既能体现电热耦合特征又不会让MILP求解时间失控。数据方面风电出力曲线用某风电场实际出力数据做归一化处理热负荷曲线按季节典型日构造。历史样本取120天每天24个时段这样原始数据矩阵是24×120。5.2 对比方法设计一套完整的对比实验应该包括以下四类方法方法不确定性处理方式特点确定性模型取预测值不考虑不确定性成本最低但不可行风险最高传统随机规划用经验分布p₀K个场景依赖精确概率假设传统鲁棒优化box不确定集最坏情况成本最高过于保守本文DRO模糊集内最坏分布折中方案用这四类方法分别求解同一系统对比总成本、弃风率、热负荷削减率等指标就能直观看到DRO在稳健性和经济性之间的平衡效果。关键的一点是要增加样本外测试环节。具体做法是用前120天的数据构造模糊集并求解再用后30天未参与训练的数据作为真实场景进行回代检验统计约束违反率和实际运行成本分布。样本外测试能有效防止过拟合到历史数据的假象这是审稿人非常关注的一点。5.3 结果分析的角度与图表我每次做结果分析必看三个角度第一成本-稳健性帕累托曲线。固定模糊集半径从0变化到0.2画出总成本随半径的变化曲线。半径越大成本越高曲线越陡说明系统对分布偏差越敏感。这套曲线能直接回答模糊集取多大合适这个实际问题。第二最坏分布的结构分析。把CG迭代最终收敛到的最坏分布p*与经验分布p₀做对比观察哪些场景的概率被上调了。通常概率上调的是风电出力偏低、热负荷偏高的恶劣场景这符合物理直觉也能验证模型确实识别出了风险。第三迭代收敛曲线。画出上下界随迭代次数的变化CG的收敛曲线通常是前几轮快速收敛后面逐渐平稳。如果出现锯齿状震荡要回去检查切割平面是不是加错了。6. 工程化落地中的坑与经验最后这部分把我实际跑这个项目时踩过的坑和积累的经验整理出来按重要性排序。6.1 求解器选择与数值稳定性主问题是MILP我用的是GurobiCG迭代配合YALMIP接口非常顺。如果暂时没有Gurobi的许可用免费的SCIP或CBC也能跑但求解速度会差不少尤其当储热罐的充放热二进制变量超过100个时差距很明显。数值稳定性方面要特别注意量纲。电功率是MW级成本是万元级如果变量尺度相差过大求解器的数值容差会引起奇怪的不收敛现象。我习惯把成本统一折算成万元、功率统一折算成MW并在YALMIP里设置求解器的数值容差参数。另外第二阶段LP里如果出现退化情况即多个最优解对偶变量的取值可能不稳定导致CG切割平面质量下降。解决办法是在目标函数里加一个极小的正则项比如0.0001倍的变量平方和实测能有效稳定对偶变量。6.2 场景数目与模糊集半径的调参经验场景数K和模糊集半径θ是两个最核心的调参对象。从我的大量测试看二者存在一定的替代关系K增大时经验分布p₀更接近真实分布所需的最小θ可以相应减小反之K较小时需要更大的θ来覆盖分布估计误差。给一个粗略的参考范围K30时θ₁取0.03到0.08通常有较好的样本外表现θ₁超过0.15后模型就趋近于传统鲁棒优化了失去了分布鲁棒的灵活性。调参时可以做一个二维扫描K×θ画出样本外成本的等高线图找到谷底区域这是最稳的方法。6.3 收敛判据与迭代加速技巧CG的收敛判据一般看相对间隙(UB-LB)/UB 0.01。但注意如果模糊集半径很小子问题F(x)随x变化本来就平缓上下界差距很小容易提前收敛到局部早停反之如果半径过大子问题波动大可能需要额外几轮迭代。一个很实用的加速技巧是热启动。第一轮求解主问题时把上一轮迭代得到的最优x作为初始可行解传给求解器能明显减少MILP求解时间。YALMIP里可以用assign和sdpsettings(usex0, 1)实现。提示如果迭代过程中出现主问题无解不要急着改代码先检查第二阶段约束定义是否有笔误尤其是h_k和T_k的维度是否匹配。这类问题比算法问题频繁得多。6.4 从能跑到可信的最后一公里代码能跑出结果之后一定要做两组验证。第一组是把模糊集半径设为0这时DRO应当退化为传统的随机规划结果应该与直接用p₀求解的随机规划完全一致。如果对不上说明切割平面或主问题约束有错。第二组是把K个场景中的某几个场景的概率设为固定值并收紧模糊集结果应当与已知的解析解或商业软件结果一致。这两组验证通过代码才算是真正可信。我在实际项目中发现很多复现论文代码的人摔倒在半径设0退化验证这一步。原因大多是主问题里加入的最坏场景切割没有正确覆盖所有场景或者概率变量p没有参与到切割构造中。排查方法很简单在CG迭代结束前把当前x代入原始DRO问题直接计算目标值跟主问题目标值对比如果偏差超过阈值说明切割缺失。从我个人经验来说分布鲁棒优化的Matlab实现难点从来不在某个数学步骤本身而在把各个模块正确组装并验证。只要守住先退化验证、再算例分析、最后样本外测试这条流程你复现的代码质量就超过大多数论文附带的源码了。现在这篇文章里给出的框架已经足够支撑你在这个方向上独立做出一套可发表级别的实验结果。