数据驱动的两阶段分布鲁棒优化在电热综合能源系统调度中的应用

这两年做综合能源系统优化，我几乎把所有的休假时间都耗在了“电热耦合 + 不确定性”这两个词上。如果你也正在做或者准备做电热综合能源系统的调度研究，那你大概率会遇到同一个问题：风电和热负荷都有随机性，传统的随机规划需要精确概率分布，现实中拿不到；纯鲁棒优化又保守得让你不敢把结果写进论文和报告里。我自己的做法，也是这次要分享的核心，就是“数据驱动的两阶段分布鲁棒优化，用1-范数和∞-范数去约束分布模糊集”，配套Matlab代码实现。这其实就是给不确定性建模套上一个由历史数据标定的“分布区间”，然后在最坏概率分布下做两阶段调度决策。

这篇内容适合正在做综合能源系统优化方向的研究生、工程师，以及所有被不确定性和调度协调问题折磨过的同行。我会把从模型推导到Matlab代码落地的完整链路讲清楚，尤其是CCG算法的迭代逻辑、不同范数约束的决策差异，还有那些论文里不会写但你调试代码时一定会踩的坑。

1. 问题拆解：电热系统为什么需要“两阶段分布鲁棒”

1.1 “以热定电”的困局：电热耦合是核心矛盾

先聊一个最基本的常识。在北方的冬季供暖场景里，热电联产机组（CHP）是主角，它的电出力和热出力之间有一条强耦合的运行区间约束。简单说就是：热负荷上来了，机组就得使劲产热，而产热的同时电出力也被“绑定”在一个比较高的水平上，这就是常说的“以热定电”。可问题是，夜间风电大发、热负荷也高，CHP又降不下来电出力，风电就被压着没法上网，形成弃风。

我刚做这个方向时，脑子里第一个想法就是：能不能用电锅炉、热泵、储热罐去解耦？理论上当然可以，但设备一多，调度决策的维度就上来了。电锅炉可以把多余的风电转化为热能替代一部分CHP的热出力，从而为CHP腾出降电出力的空间。可这又依赖预测数据：风电预测偏差一上来，调度方案就不可行。所以，不确定性的处理就变成绕不过去的一环。

1.2 不确定性处理：为什么最终选了分布鲁棒

在接触分布鲁棒之前，我其实先试过随机规划和鲁棒优化。随机规划要求你给出不确定变量的精确概率分布，然后用场景期望去衡量目标函数。问题是，在实际工程中，历史样本再多，也只能估计出一个近似分布，用近似的分布去优化，结果就可能“差之毫厘，谬以千里”，尤其在分布尾部风险上特别敏感。

鲁棒优化则走向另一个极端，它假设不确定变量落在某个确定集合内，并针对最坏情况做决策。好处是绝对稳健，坏处是调度成本高到离谱，甚至在经济性上根本不可接受。

分布鲁棒优化恰好卡在两者中间。它的核心思路是：不再假设不确定变量服从某个精确分布，而是构造一个包含“所有可能概率分布”的模糊集，然后在这个模糊集内寻找最坏情况下的期望成本，并对此做决策。而这个模糊集的构造，正是通过历史数据得到的。换句话说，数据越多、越可靠，模糊集就越小，决策就越贴近真实最优；数据越少，模糊集越大，决策就越保守。这就是“数据驱动”四个字的真正含义。

我最终选择1-范数和∞-范数约束的原因也很直接：这两种范数构造的模糊集是线性的，可以用线性规划或混合整数线性规划求解，配上列与约束生成算法（C&CG），在Matlab里实现非常顺滑。相比之下，如果用Wasserstein距离去构造模糊集，虽然理论上更漂亮，但计算复杂度会明显上升，对算力有限的场景不够友好。

2. 模型构建：两阶段框架与模糊集设计

2.1 两阶段优化的完整框架

两阶段分布鲁棒优化的名字听起来唬人，但拆开就很好理解。第一阶段对应的是“日前决策”，也就是在不确定性还没发生的时候，你需要提前定下来的变量。在电热系统里，就是CHP机组的启停状态、储热罐的充放热计划、以及一些设备的投运状态。这些变量有一个特点：一旦定了，短时间内改不了，或者改动成本非常高。

第二阶段对应的是“实时调整”，也就是不确定性实现之后，你可以做的修正动作。比如调节CHP的实际电出力、热出力，增加或减少弃风量，必要时切一些负荷，调整电锅炉的功率等。这些变量可以在看到风电和热负荷的真实值后再决定，代价是偏离日前计划会产生调整成本，比如考核费用、煤耗变化等。

整个目标函数也就顺理成章拆成两部分：

目标函数 = 第一阶段成本（日前决策）+ 第二阶段最坏期望成本（实时调整）

这里的“最坏期望”怎么理解？设想你手里有一堆历史数据和由这些数据构造出来的、包含多个可能概率分布的模糊集。你在这个模糊集里找到一个概率分布，使得第二阶段调整成本在所有分布中是最高的，而你的优化目标就是让这个“最高成本”尽可能小。这就是典型的min-max-min结构，也是分布鲁棒优化的核心数学形式。

2.2 1-范数和∞-范数约束的模糊集到底在约束什么

现在到了这个方向最核心、也最容易被误解的地方：1-范数和∞-范数约束的模糊集，到底在约束什么？

先说基本设定。假设你有K个历史场景，每个场景对应一组风电出力和热负荷的样本。初始时，每个场景的概率可以取等概率，也就是p0_k = 1/K。但是历史数据只是真实分布的抽样，真实概率分布和这个经验分布之间一定存在偏差。那怎么量化这个偏差呢？我们用1-范数和∞-范数来度量。

• 1-范数约束：所有K个场景的概率偏差绝对值之和不超过θ1，也就是 (\sum_{k=1}^{K} |p_k - p0_k| \le θ1)。
• ∞-范数约束：所有K个场景中，单个场景的概率偏差绝对值最大值不超过θ∞，也就是 (\max_{1 \le k \le K} |p_k - p0_k| \le θ∞)。

换个方式理解。1-范数约束的是“整体水平”，它要求所有场景概率偏差加起来不能太大，保证模糊集整体不会偏离经验分布太远。∞-范数约束的是“单点水平”，它盯住某一个概率偏差最大的场景，防止某一个场景的概率被高估或低估得离谱。

两者之间还有一个重要的等价比率关系：(\theta_\infty \le \theta_1 \le K \theta_\infty)。这意味着，当你单独用1-范数约束时，有可能出现某个单场景概率偏差特别大的情况；而你单独用∞-范数约束时，所有场景加起来的总偏差可能超出你预期。所以工程实践里最常见的做法是同时使用这两种约束，让模糊集既在整体上受限，又在单点上受限。

θ1和θ∞怎么确定？这不是拍脑袋定的，而是基于历史样本数N和置信度β推导出来的。常见的公式是：

[ \theta_1 = \frac{K}{2N} \ln\frac{2K}{1-\beta},\quad \theta_\infty = \frac{1}{2N} \ln\frac{2K}{1-\beta} ]

从公式里可以明显看出一个趋势：样本数N越大，θ越小，模糊集越“瘦”，结果越贴合真实；置信度β越高，θ越大，模糊集越“胖”，决策越保守。这个性质在后面做敏感性分析时会反复遇到。

2.3 完整模型与约束写法

我直接给出一套我实际用过的、可复现的简化模型写法，方便大家和后面的代码对上号。注意这是一个面向调度的经济性模型，不涉及网络潮流细节，适合先跑通算法再做扩展。

第一阶段决策变量：(x) 表示CHP启停、储热罐状态等离散和连续变量的集合。目标函数包括启停成本、日前计划出力成本。

第二阶段决策变量：(y) 表示在给定任意一个场景ξ（风电、热负荷）后的调整量，包括机组出力调整、弃风量、切负荷量、储热罐调整量等。

模型整体形式可以写成：

[ \min_{x} \left[ c^T x + \max_{p \in \Psi} \sum_{k=1}^{K} p_k \min_{y_k} d^T y_k \right] ]

其中 (\Psi) 就是由1-范数和∞-范数约束构造的模糊集。

约束条件如下：

• 电功率平衡约束：第一阶段电出力 + 风电 + 第二阶段调整量 = 电负荷。
• 热功率平衡约束：CHP热出力 + 电锅炉/热泵产热 + 储热罐放热 + 调整量 = 热负荷。
• CHP可行域约束：电出力和热出力必须落在多边形的运行区间内，这是“以热定电”的数学化身。
• 储热罐约束：储热量的上下限、充放热速率限制、以及日始日终储热量一致约束。
• 弃风与切负荷惩罚约束：在第二阶段允许一定量的弃风，但必须付出惩罚成本；同理切负荷也有限额和惩罚。

这套模型的完整程度，已经足够支撑一篇高质量论文的算例，也足以让你理论上理解电热协调调度的所有关键环节。接下来最关键的，就是怎么把它解出来。

3. 求解算法与Matlab代码实现

3.1 为什么不能直接用一个优化求解器搞定

很多第一次接触两阶段分布鲁棒的同学，拿到模型后第一反应是把整个模型扔给Yalmip，让它一把梭哈。结果往往是要么内存爆炸，要么求解时间长得让人绝望。原因很简单：这个问题是min-max-min三层结构，不是一个标准线性规划可以表达的。

我用的求解方法是列与约束生成算法（C&CG），这是目前在两阶段分布鲁棒和两阶段鲁棒优化里最主流的方法。它的思想本质上是“动态添加约束”，一边迭代一边把关键约束补进主问题里。

你可以把它想象成谈合同：主问题先提出一个相对宽松的计划，子问题去“找茬”，找到最坏场景下的违约点，再把新条款写进合同里，继续谈。反复几轮之后，合同条款足够全面，两边就达成一致了。

3.2 C&CG算法的完整迭代流程

C&CG算法每一步做什么，我拆开来讲。

第一步，初始化。下界LB设为负无穷，上界UB设为正无穷。选择一个初始的场景，通常是历史样本里的第一个或均值场景，作为子问题的初始不确定性参数。

第二步，求解主问题。主问题是一个包含现有所有割约束的min问题，求解后得到第一阶段决策x和一个下界更新值LB。这里的主问题本质上还是混合整数线性规划，因为第一阶段有启停状态这种0-1变量。

第三步，固定x，求解子问题。子问题是max-min结构：先内层求min（给定一个概率分布和场景下的第二阶段最优成本），外层求max（找最坏概率分布）。实际操作中，我会把内层min问题先对偶成max问题，这样内外层都是max，合并成一个单层max问题，然后用求解器直接解。

第四步，更新上界。子问题求得的最坏期望成本加上第一阶段成本，就是当前决策下的总成本，这个值可以作为UB。

第五步，收敛判断。如果UB和LB的相对间隙小于设定阈值，就停下来，输出当前的x和y作为最优解。否则，把子问题找到的最坏场景对应的第二阶段约束写成新的割约束，加入主问题，然后回到第二步。

需要特别注意的是，CCG和Benders分解的核心区别：Benders是给原问题加关于目标值的割平面，而CCG是直接把新的场景约束加进去，相当于“跳过了”中间层，所以收敛速度明显更快。我实测下来，CCG通常在10到50次迭代内就能收敛到1e-3的相对间隙，而Benders可能要上百次。

3.3 Matlab代码实现框架与关键代码

数学层面有了，接下来就是真正动手写代码。我的环境是Matlab + Yalmip + 商业求解器，这套组合在学术界和工程界都很常见，优点是建模方便，迭代循环里动态加约束也好实现。

先给出代码的整体框架：

%% 1. 参数定义 % 系统参数：节点数、CHP参数、储热罐参数、历史场景数据等 load('data_history.mat'); % 历史风电/热负荷场景，K行N列 K = size(data_history, 1); % 场景数 %% 2. 模糊集参数计算 beta = 0.95; % 置信度 theta_inf = 1 / (2*N) * log(2*K / (1-beta)); theta_1 = K / (2*N) * log(2*K / (1-beta)); %% 3. CCG主循环 LB = -inf; UB = inf; iter = 1; while abs((UB - LB)/LB) > 1e-3 && iter < 100 % 求解主问题：mixed-integer LP ops = sdpsettings('solver', 'cplex', 'verbose', 0); optimize([cons_main, cuts(1:iter-1)], obj_main, ops); LB = value(obj_main); x_val = value(x); % 固定x，求解子问题：max-min -> 对偶后单层max [sub_obj, worst_profile] = solve_subproblem(x_val, theta_inf, theta_1); UB = c' * x_val + sub_obj; % 生成新割约束 if abs((UB - LB)/LB) > 1e-3 cuts(iter) = generate_cut(worst_profile); iter = iter + 1; end end

主问题的建模，简化之后有下面几个核心块。Yalmip里定义变量、加约束、写目标，然后调用求解器。这里面有一个特别容易出错的地方：割约束是动态增长的，所以在循环外要预留一个空的cut容器，每次迭代往里塞新的约束。

% 主问题建模示意 x = binvar(n_units, 1); % CHP启停 p_chp = sdpvar(n_units, 1); % CHP电出力 h_chp = sdpvar(n_units, 1); % CHP热出力 % ... 其他变量 obj_main = cost_start(x) + cost_fuel(p_chp) + eta; % eta是第二阶段成本的下界变量 cons_main = [power_balance(p_chp, wind, load), ... heat_balance(h_chp, boiler, storage), ... chp_feasible_region(p_chp, h_chp, x), ...];

子问题，也就是第二阶段，是最麻烦的部分。我直接把max-min合并成对偶之后的max问题。这里的关键一步是把第二阶段的min问题写成对偶形式。对偶之后，内层min被消掉，变成外层max对偶变量加目标函数，同时要满足对偶可行约束。

function [obj, worst_profile] = solve_subproblem(x_val, theta_inf, theta_1) % 定义第二阶段变量 y 和对偶变量 % 核心：maximize sum(lambda_k * (h - H*x_val)) + ... % 约束：对偶可行、范数约束 % 返回：最坏场景profile和子问题目标值 end

有一点必须提醒：子问题的场景数K直接决定了对偶问题的变量和约束规模。场景数太大，比如超过1000个，求解时间会明显拉长。我自己的做法是先做场景削减（比如用K-means聚类把历史场景聚成50～100个有代表性的），再进模糊集，效果很好。

3.4 求解器配置与仿真环境

代码本身要依赖求解器。我的建议是至少配置一个能够解决混合整数线性规划的求解器。单纯用Matlab内置的linprog或intlinprog也可以做小规模算例，但一旦系统节点数到几十个、场景数过百，内置求解器的效率就很难看了。

在Yalmip里，选择求解器就一行命令：

sdpsettings('solver', 'cplex', 'verbose', 2);

其中solver可以替换成你本机安装的求解器名称。跑之前务必用Yalmip自带的yalmiptest命令验证一下求解器是否配置成功，这是最省时间的一步，能避免浪费两小时在环境配置上。

4. 算例结果分析：不同范数约束之下的决策差异

4.1 算例设置

我用的是一个经典的6节点电力/6节点热力系统（指代某公开的测试系统），系统里包含两台CHP机组、一台电锅炉、一个储热罐和一个风电场。历史数据用随机性建模的方法生成300个场景，覆盖风电出力和热负荷的联合分布。基准模型的置信度设为0.9，样本数取300，场景削减后选50个代表性场景。

为了验证模糊集的有效性，我设计了三组对比：

每一组都跑CCG算法，记录总成本、迭代次数和最终的弃风量。

4.2 核心结果与解读

先看总成本。随机规划的总成本最低，这个没有悬念，因为它默认你掌握了真实分布，不存在“怕出事”的额外代价。1-范数约束的分布鲁棒模型，总成本比随机规划高大约4%到6%，这部分多出来的成本可以看作“买保险”的费用。加上∞-范数约束之后，总成本再上浮1%到2%，但换来的是最坏场景下的成本明显下降。

再论弃风量，这是电热系统最直观的指标。随机规划的最优解里，弃风量在多数场景下都是0，但是一旦风电预测偏差朝不利方向偏移，弃风量飙升。分布鲁棒模型的弃风量在各个场景下的分布更平滑，尤其是加入∞-范数后，极端场景下的弃风被明显抑制了。折线图上看，随机规划的“最高点”很高，M3的“最高点”则被压扁。

这里要特别说明一下1-范数和∞-范数在实际表现上的差异。1-范数约束对“总体的概率偏差”敏感，它会保证大部分场景的概率分布不至于太离谱。但是，它可能允许某一个特别关键的场景概率被高估。而∞-范数约束恰好补上这个漏洞——它把单个场景的概率偏差封死在一个小范围内。所以同时用两者，相当于“既看整体，又抓重点”。

4.3 样本量与置信度的敏感性分析

这个部分我建议每个做这类模型的人都跑一遍。固定其他参数不变，只改变历史样本数N从50到500，观察目标函数值和θ的变化。结果是：目标成本随N增大而下降，并且下降速率在N小于200时很快，之后趋于平稳。原因就是样本足够多之后，经验分布已经非常接近真实分布，模糊集的收紧空间有限了。

置信度的敏感性方向正好相反。β从0.8升到0.99，目标成本持续上升，因为更高的置信度意味着更大的θ，分布“不确定的范围”变大，系统必须为更极端的情况留出备用。在写论文时，这个指标可以用来说明“决策者的风险规避程度和经济性之间的权衡”。

还有一个所有做这类算法的人都会关心的指标：迭代次数。我实测了几组场景下的CCG收敛情况：

可以看到，场景数翻倍，迭代次数只是微增，但子问题求解时间几乎成线性增长。所以，如果你追求快，就尽量通过场景削减把K控制在100以内。

5. 踩坑记录与实操建议

这个部分是重点中的重点。以下五个坑是我在实际调代码过程中真正踩到过的，而且每一个都让我多花了好几天时间。

5.1 最坑的五个实操问题

第一个坑是初始概率p0的设置。很多人一上来就把p0设为等概率1/K，觉得这是理所当然的。问题在于，模糊集里的p被约束在p0附近，如果p0本身和真实分布差得很远，整个模型的“数据驱动”效果就会大打折扣。建议的做法是：先对历史场景做聚类，按聚类簇的占比作为初始概率，再在这个基础上构建范数模糊集。这样模糊集的初始位置就更准了。

第二个坑是子问题对偶时的大M参数设置。把第二阶段min问题对偶成max之后，需要把一些线性化的双线性项处理掉，通常会引入大M。大M设太大，求解器会出现数值问题，对偶变量失控；设太小，又可能把真实可行域截断。我自己的经验是让大M取系统最大功率的10倍，然后根据结果微调。

第三个坑是Yalmip里双线性项的处理。如果你直接在Yalmip里写p_k乘以y_k，它会把问题识别成非线性非凸问题，求解器直接罢工。解决办法是，要么在数学层面先做一个变量替换处理掉双线性项。要么在代码层面把子问题拆成两层循环，用穷举场景和逐场景求解的方式绕过这个坑。我在最初调试时走了弯路，最后是把子问题变成对偶问题后，所有双线性项都自然消失了。

第四个坑是割约束的更新逻辑。CCG的割约束生成时，必须和主问题里已有的变量建立关联，不是简单地把一组数塞进去。具体来说，新增的割约束是“第二阶段可行性和最优性的KKT条件”在某个具体场景下的投影，它引用的是主问题里的x变量。如果写错了关联关系，CCG会陷入死循环——LB和UB几千次迭代都不收敛。

第五个坑是求解器选择与数值容忍度。整数变量和连续变量混在一起时，求解器的MIP容忍度设置会影响迭代次数。默认的容忍度通常是1e-4，但在CCG里，我建议把它放宽到1e-3，收敛速度会明显加快，而结果精度几乎没有下降。因为CCG本身也会在间隙小于1e-3时停止，过高的内部精度只会白白增加求解时间。

5.2 参数调优与收敛性建议

如果你跑出来的迭代曲线迟迟不收敛，我建议按下面的顺序排查。

第一步，检查θ的值是否过大。θ如果比0.5还大，说明模糊集已经“膨胀”到失去了限制意义，这时候先降低置信度，或者增加历史样本数。

第二步，缩小第二阶段模型的复杂度。第二阶段变量越多、约束越复杂，对偶问题也就越难解。可以先把储热罐的候选状态从24小时简化为3个典型时段，验证CCG收敛后再恢复完整模型。

第三步，观察每轮迭代的最坏场景是否在反复横跳。如果新生成的割约束对应的场景总是在两个极端之间来回切换，说明子问题在数值上有振荡。这时可以在子问题里增加一个正则化项，比如给弃风量加一个极小的二次惩罚系数。

最后，我个人还有一个习惯：每轮CCG迭代之后，把主问题解出的x存下来画图。如果发现两个相邻迭代之间x的变化在震荡，别怀疑，肯定是割约束写错了，回去看第五个坑。

6. 写在最后的一些心得

模型和方法说完了，再分享一点我个人的真实感受。两阶段分布鲁棒这个方向，入门难在数学推导，进阶难在代码实现，但真正的瓶颈其实是“把模型讲清楚”。我自己最初读到这类论文时，最大的感受是：公式密密麻麻，但读者很难抓住“1-范数到底带来了什么变化”。所以我在自己的研究和代码里，特别注重做对比实验——用同一个算例，分别跑随机规划、鲁棒优化、1-范数、∞-范数和混合范数，把成本和弃风情况摆在桌面上。没有这些对比，分布鲁棒的价值是说不清的。

另外，我知道很多人会从论文里直接抄一个算例来跑代码，但我的建议是：先在自己的小系统上，比如6节点加6节点的热网这个规模上把CCG跑通，再往更大、更复杂的系统迁移。因为CCG算法本身的收敛性和数值稳定性，足够你在小规模上摸透所有的坑。等算法逻辑完全清楚了，再换IEEE标准算例或者实际工程数据，才不会一头扎进调试地狱。

最后再分享一个可以扩展的方向：这套两阶段分布鲁棒框架不只能用在电热系统上，只要你能把不确定性写成概率分布模糊集，它几乎可以平移到任何含不确定性的调度问题上，不管是电气综合能源、交通供电系统还是多微网协同。而且用1-范数和∞-范数的好处是，整个模型始终保持在混合整数线性规划的框架内，扩展起来逻辑非常统一。往后的研究里，你大可以把多阶段扩展、碳交易机制、需求响应全部嵌进来，框架仍然成立。

如果看完这些，你对CCG迭代或模糊集参数的推导还有疑问，欢迎留言交流。我自己当初也是花了很多个不眠夜才把这些逻辑彻底理顺的，所以特别愿意看到后来者能少走一点弯路。