Matlab实现两阶段分布鲁棒优化：电热综合能源系统调度实战

近年来在做电热综合能源系统调度研究时，我明显感觉到一个趋势：纯物理建模的路子越来越走不动了。原因倒不复杂——热网动态特性、建筑热惯性、可再生出力波动这些环节，实际运行数据里藏着大量模型难以精确刻画的信息，尤其在样本量不足的小场景下，传统数据驱动模型容易过拟合，而纯物理模型又泛化不足。这篇要分享的Matlab实现，用的是数据驱动的两阶段分布鲁棒优化思路，通过1-范数和∞-范数约束来构造不确定性集合，把电热综合能源系统的日前调度问题转化为一个可求解的数学规划模型。

这套方法的核心价值在于：它不假设不确定参数服从某个固定分布，而是从历史数据中提取统计信息，构造一个包含真实分布在内的分布模糊集。1-范数约束保证了模糊集不会过松，∞-范数约束保证了不会过紧，两阶段框架则把“现在决定”的机组组合问题和“等不确定性实现后再调整”的经济调度问题分开处理。对于正在做综合能源系统优化、分布鲁棒优化或者Matlab建模仿真的同学，这篇内容应该能帮你省下不少绕弯子的时间。

我会从问题建模、不确定性集合构造、两阶段模型的求解思路、Matlab代码实现细节到算例分析，一步步拆开来讲，顺便把我自己在实际调试中踩过的坑一并列出来。

1. 为什么偏偏是“两阶段分布鲁棒”：问题动机与建模思路

1.1 电热综合能源系统的调度难点在哪里

电热综合能源系统（Integrated Electricity and Heating System, IEHS）的核心矛盾在于电和热两种异质能源的强耦合。热电联产机组（CHP）是这种耦合的典型代表——它同时产电和产热，但电出力与热出力之间存在可行域约束，不能像独立电厂那样随意调节。再加上热网本身的传输延迟和储热罐的缓冲作用，系统的调度决策就变成一个跨时间、跨介质、含不确定性的复杂优化问题。

传统做法是确定性优化：给定风电预测出力曲线和负荷预测曲线，求解一个最小化运行成本的机组组合问题。问题在于，风电的实际出力与预测值偏差在冬季供暖期经常能达到20%-30%，预测误差一旦变大，按照确定性方案执行的调度结果轻则经济性变差，重则导致切负荷或弃风。

鲁棒优化是另一种思路，它把不确定参数放在一个确定的集合里，然后求解最坏情况下的最优决策。但传统鲁棒优化的结果通常过于保守——它要求所有不确定参数同时取到最坏值，这在现实中几乎不会发生。

1.2 数据驱动分布鲁棒与“固定分布”思路的本质区别

分布鲁棒优化（Distributionally Robust Optimization, DRO）的思路介于两者之间：它不假定不确定参数服从某个精确的概率分布，而是利用历史数据构造一个包含真实分布在内的分布模糊集（Ambiguity Set），然后求解在最坏分布下的期望成本最小化。

这样做的好处很直观：当历史数据充足时，模糊集收缩，模型逼近随机规划；当数据不足时，模糊集扩大，模型趋向鲁棒优化。这种“由数据决定保守程度”的特性，让DRO特别适合风电出力这类分布不确定性强、历史样本有限的对象。

1.3 两阶段框架在电热调度里到底扮演什么角色

两阶段分布鲁棒模型对应到实际业务场景，可以这样理解：

• 第一阶段（Here-and-Now）：在不确定性实现之前，决定机组的启停状态、CHP的产热计划、储热罐的充放热计划等“慢变量”。这些决策需要提前确定，一旦执行难以快速更改。
• 第二阶段（Wait-and-See）：当风电实际出力、负荷实测值等不确定参数揭示后，根据第一阶段决策，调整各设备的实际出力水平，目标是在满足约束的前提下把调整成本降到最低。

两阶段的好处是决策层次清晰：第一阶段保证系统可行，第二阶段保证经济最优。从数学上看，这个模型是一个min-max-min结构的优化问题，外层最小化总成本，中间层对最坏分布求期望，内层对每个场景求二阶段调整成本，求解难度比单阶段模型上了一个台阶，需要通过强对偶理论或KKT条件将其转化为单层优化问题。

2. 1-范数和∞-范数约束到底在约束什么：不确定性集合的数学构造

2.1 从历史数据中构造分布模糊集的步骤

假设我们有S个风电出力的历史样本，每个样本对应一个离散场景ξ_s。数据驱动DRO的做法是：以这S个样本的经验分布为基准，构造一个“真实的分布”可能落入的集合。

在实现时，我们用概率向量p = [p_1, p_2, ..., p_S]来表示每个场景的概率权重，经验分布对应p_0 = [1/S, 1/S, ..., 1/S]。模糊集就定义为：与经验分布的距离不超过某个阈值θ的所有分布的集合。

这个“距离”的度量方式直接决定了模型的保守程度和求解难度。我在代码里用的是两种经典范数约束的组合——1-范数和∞-范数。

2.2 1-范数和∞-范数约束的数学表达与保守性对比

1-范数约束的表达式为：

||p - p_0||_1 ≤ θ_1

即 Σ|p_s - 1/S| ≤ θ_1。这个约束限制的是所有场景概率偏差的绝对值之和，它控制的是分布的整体偏移程度。

∞-范数约束的表达式为：

||p - p_0||_∞ ≤ θ_∞

即 max|p_s - 1/S| ≤ θ_∞。它限制的是单个场景概率偏差的上限，控制的是局部偏差。

代到代码实现里的完整约束是：

Σ|p_s - 1/S| ≤ θ_1 max|p_s - 1/S| ≤ θ_∞ θ_1, θ_∞ 的值由历史数据量和置信度决定

两个范数配合使用的直观理解：

• 1-范数约束控制“整体漂移”。如果只用∞-范数，可能会出现每个场景的概率都小幅偏移、累积起来偏差很大的情况。
• ∞-范数约束控制“单点异常”。如果只用1-范数，可能会出现个别场景概率被大幅高估或低估的情况。

把两者结合，等于给概率分布同时上了“总量限制”和“分量限制”，模糊集在两个方向上都被锁住。从实验结果看，组合约束的保守性介于单一∞-范数和单一1-范数之间，但稳定性最好。

2.3 置信度与模糊集半径的量化关系

这里有一个很多人容易忽略的细节：θ_1和θ_∞怎么取值？取大了模型过于保守，取小了模糊集可能没有覆盖真实分布，模型失效。

从统计学的角度，如果置信度为α，历史样本数为S，θ_1和θ_∞的取值可以按下式估算：

θ_1 = (2/S)·ln(2K/(1-α))^(1/2) θ_∞ = (1/2S)·ln(2K/(1-α))^(1/2)

其中K为不确定参数的总数。这两个公式的思路是：样本越多，数据越能说明问题，模糊集半径应该越小；置信度越高，越需要把模糊集做大，保证真实分布大概率被包含在内。

实际代码中，我通常的做法是给定一个α，然后通过设置不同的θ组合来观察调度成本和鲁棒性的变化，绘制帕累托前沿，再结合工程实际确定一组折中值。

3. 两阶段三层结构的数学模型：目标函数与约束条件

3.1 第一阶段决策变量与运行约束

在Matlab代码里，我把第一阶段的决策变量定义为：

第一阶段的约束包括机组最小启停时间约束、CHP电热可行域约束、储热罐容量约束和热平衡约束。这些约束的共同特点是：它们只依赖第一阶段的决策，不涉及不确定参数的实现，因此在优化开始前就可以确定。

3.2 第二阶段决策变量与实时调整约束

到了第二阶段，风电出力ξ_s的真实值已经“揭晓”。此时，系统需要在第一阶段决策的基础上，对常规机组的电出力、CHP的电出力、电锅炉的用电功率等变量进行再调整。

第二阶段的决策变量主要包括：

• ΔP_i：常规机组i的出力调整量
• ΔP_CHP：CHP电出力调整量
• P_EB：电锅炉的实际用电功率
• P_curt：弃风量

第二阶段的约束包含功率平衡约束、机组爬坡约束和出力上下限约束。注意这里的约束都带有场景下标s，表示每个场景下都需要满足。

3.3 min-max-min结构的目标函数与对偶转换思路

两阶段DRO的完整目标函数形式为：

min {第一阶段成本 + max_{p∈Ω} Σ_s p_s · Q(x, ξ_s)}

其中Q(x, ξ_s)是第二阶段在场景s下的最优调整成本，p是场景概率分布，属于模糊集Ω。

这个三层结构没法直接丢给求解器。标准处理步骤是：

1. 把内层的max问题拆出来。由于内层是关于概率p的线性问题，而模糊集由1-范数和∞-范数约束定义，恰好是线性约束，因此内层max问题可以等价变换为一个包含对偶变量的优化问题。
2. 把第二阶段的min Q(x, ξ_s)用其对偶问题替换，与外层min合并。
3. 最终整个模型被转换为一个单层的混合整数线性规划（MILP），可以直接调用Cplex或Gurobi求解。

这里有个非常容易踩坑的地方：第二阶段问题的对偶转换要求原问题是线性的，如果加入了非线性的网损约束或非凸的可行域，整个模型就会变成难以求解的MINLP。因此在实际建模时，需要合理线性化电网网损和热网动态。

4. Matlab代码实现：从模型到可运行代码的完整路径

4.1 代码的整体架构与文件组织

我的实现用Matlab + YALMIP工具箱 + Cplex求解器。工程代码按以下结构组织：

IEHS_DRO/ ├── main.m % 主程序入口 ├── data/ │ ├── wind_data.mat % 风电历史数据 │ ├── load_data.mat % 电/热负荷数据 │ └── system_params.m % 系统参数配置 ├── model/ │ ├── build_1st_stage.m % 第一阶段约束构建 │ ├── build_2nd_stage.m % 第二阶段约束构建 │ ├── build_ambiguity.m % 模糊集约束构建 │ └── build_obj.m % 目标函数组装 ├── solver/ │ ├── solve_iehs.m % 求解主逻辑 │ └── dual_transform.m % 对偶转换函数 └── result/ └── plot_results.m % 结果可视化

4.2 YALMIP建模核心步骤与代码片段

第一步：定义变量。注意第一阶段二进制变量要用binvar，第二阶段连续变量用sdpvar，场景变量要带场景维度。

% 第一阶段变量 u = binvar(n_gen, T, 'full'); % 机组启停 x = binvar(n_gen, T, 'full'); % 开机动作 H_CHP = sdpvar(1, T, 'full'); % CHP热出力计划 S_TES = sdpvar(1, T+1, 'full'); % 储热状态 % 第二阶段变量（每个场景一组） P_adj = sdpvar(n_gen, T, S, 'full'); % 机组出力调整量 P_EB = sdpvar(1, T, S, 'full'); % 电锅炉功率 P_curt = sdpvar(1, T, S, 'full'); % 弃风量 % 概率变量 p = sdpvar(S, 1, 'full'); % 场景概率

第二步：构建模糊集约束。这里用辅助变量把绝对值约束线性化。

% 概率和为1 F = [sum(p) == 1, p >= 0]; % 1-范数约束：sum(|p - p0|) <= theta_1 d1 = sdpvar(S, 1, 'full'); p0 = ones(S, 1) / S; F = [F, -d1 <= p - p0 <= d1, sum(d1) <= theta_1]; % ∞-范数约束：max(|p - p0|) <= theta_inf dinf = sdpvar(1, 1, 'full'); F = [F, -dinf <= p - p0 <= dinf, dinf <= theta_inf];

第三步：构建第二阶段成本函数。Q(x, ξ_s)包含调整费用、弃风惩罚和切负荷惩罚。

% 第二阶段调整成本 adjust_cost = sum(sum(C_adj .* P_adj, 1), 2); curtail_cost = C_curt * sum(sum(P_curt, 1), 2); Q_s = squeeze(adjust_cost + curtail_cost); % 每个场景一个成本值

第四步：最关键的一步——对偶转换。在Matlab里我们不需要手推对偶问题，而是利用YALMIP有限制的对偶能力或直接构造扩展模型。一个更稳健的做法是：先写出内层max问题的对偶，再将其并入主问题。我在代码里采用了列与约束生成（C&CG）算法，避免一次性把大规模问题直接塞给求解器。

4.3 C&CG主问题与子问题的迭代求解策略

C&CG算法的核心是“主问题-子问题”交替迭代：

• 主问题（MP）：在有限个已生成的最坏分布点下求解第一阶段决策和总成本。
• 子问题（SP）：固定第一阶段决策，求解最坏分布下的第二阶段成本，并把该分布对应的最优性割平面返回给主问题。

伪代码如下：

初始化：设定最大迭代次数，UB = Inf， LB = -Inf， k = 1 循环直到 UB - LB <= 阈值： 1. 求解主问题，得到第一阶段决策x_k和目标值obj_MP LB = max(LB, obj_MP) 2. 固定x_k，求解子问题，得到最坏分布p_k*和成本obj_SP 计算目标值UB_k = 第一阶段成本 + obj_SP UB = min(UB, UB_k) 3. 如果收敛则跳出 4. 将新生成的最坏分布p_k*作为新的场景加入主问题 5. k = k + 1

这个算法的优点是不需要显式写出对偶问题的全部细节，每次迭代只需求解当前规模下的优化问题，代码结构清爽，而且收敛速度在实际算例中表现很好，通常在10次迭代以内就能达到10^-3级别的间隔。

4.4 求解器选择与参数配置注意事项

我用的是Cplex，但在YALMIP环境下Gurobi也可以。两者在MILP求解性能上差距不大，关键是要设置好参数：

ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.mip.tolerances.mipgap = 1e-4; ops.cplex.mip.tolerances.integrality = 1e-5; ops.cplex.mip.strategy.search = 2; ops.cplex.timelimit = 7200;

这里特别提醒一句：MIP Gap的容差设置要谨慎。我一开始图省事设成1e-2，结果收敛后调度成本偏离了真实最优值将近3%。对于研究论文来说这个误差可能还能接受，但如果要做工程应用，建议至少到1e-4。

5. 算例设计与结果分析：从风电数据到调度策略

5.1 系统参数与历史数据设计

我用的算例是一个改造后的IEEE 6节点电力系统与6节点热力系统耦合的测试系统。系统包含：

• 2台常规火电机组
• 1台CHP机组（背压式）
• 1个风电场
• 1个电锅炉
• 1个储热罐
• 1个热负荷节点和1个电负荷节点

历史风电数据取某风电场冬季3个月的出力记录，采样间隔1小时，共约2160个点。考虑到DRO对样本量的敏感性，我分别测试了S = 50、100、200三个样本规模，观察模糊集半径与模型保守度的变化。

5.2 不同范数约束下的调度结果对比

先看只使用单一范数约束的结果。

可以看到，单一∞-范数的模型最保守，总成本最高。组合约束的总成本只比单用1-范数高一点点，但能有效限制单个场景的概率异常，实际调度方案更稳健。

从物理角度看，这个结果也合理：1-范数约束在控制整体分布偏移方面效率更高，而∞-范数约束补上了它对单点概率失控的短板。两者组合能在大约相同总成本下，提供更强的分布鲁棒性保证。

5.3 迭代收敛曲线与求解耗时表现

C&CG算法在S=100的场景数量下，迭代9次收敛，总耗时约186秒。每轮迭代的主问题规模虽然逐步增大，但增加的都是线性约束和连续变量，MILP部分增长有限。

在实际运行中，耗时最长的不是主问题，而是子问题中第二阶段LP的批量求解。我对这部分做了场景并行化处理，利用Matlab的parfor对场景循环并行计算，子问题求解时间从原来的约15分钟压缩到不到3分钟。如果你也在做类似研究，这一步优化非常值得投入时间。

6. 实操中踩过的坑：对偶转换与约束线性化的细节

6.1 概率变量的非负约束丢失导致的对偶模型错误

这个坑是我印象最深的。第一次做对偶转换时，我手工推导后把概率变量的非负约束p ≥ 0给丢了，导致生成的模糊集实际上允许负概率出现。症状表现为：主问题迭代到某一步后，子问题的目标值出现负的调整成本，而且越迭代越离谱。

排查半天才发现问题出在对偶转换时对偶变量符号方向搞反了。后来我长了个记性：对偶变换后一定要先做小规模数值验证。具体做法是：取一个S=3的小算例，手工列出原问题和转换后问题的KKT条件，逐项检查约束对应关系是否完整。

6.2 热网动态约束线性化过程中出现的数值病态

电热综合能源系统的热网动态方程是偏微分方程，在调度模型里通常采用节点法或有限差分法离散。问题在于，当离散时间步长取的比较小时，热网管道传输延迟会引入大量二进制变量，导致模型求解时间爆炸。

我采用的替代方案是：把热网的传输延迟近似为多阶惯性环节的离散状态空间模型，舍去热损非线性的高阶项，只保留一次线性近似。这样处理之后模型规模显著下降，而调度结果的偏差控制在2%以内，在实际工程中可以接受。

6.3 储热罐约束的时序耦合处理

储热罐的状态变量在两个时间段之间存在耦合约束：

S_TES(t+1) = S_TES(t) + H_TES_in(t)·η_ch - H_TES_out(t)/η_dis

这个约束本身是线性的，但它把全时段的决策变量连在一起，导致YALMIP建模时约束矩阵的稀疏性下降。我优化了约束的添加顺序：先把所有时段的状态转移约束一次性添入，再添加设备容量约束，最后添入耦合约束。实验证明约束顺序调整后Cplex的预求解（Presolve）效率提升约20%。

6.4 风电场景生成时消除样本间相关性的经验

最终算例里的场景不是直接使用原始风电出力数据，而是先做了PCA降维和相关性去除。原因在于，原始数据中相邻时段的风电出力高度相关，直接作为独立场景输入会高估系统的调节能力。

处理方法：对风电出力向量做主成分分析，保留前5个主成分（累计方差贡献率约92%），在降维后的主成分空间内进行K-means聚类，生成代表性场景并统计每个场景的经验概率。这样得到的场景集合在保持原始数据统计特征的同时，大大减少了场景数量，降低了模型维度。

7. 模型扩展思路：从单区域向多区域与多能互补方向推进

如果想把这套方法应用到你自己的课题中，有几个扩展方向值得考虑：

• 多区域互联：把单区域的电热耦合扩展为多个区域通过联络线和热网管道互联的情形，模糊集的构造逻辑不变，但需要额外处理区域间的耦合约束，C&CG子问题的规模会成倍增长。
• 加入碳交易机制：在目标函数中增加碳配额成本和碳交易成本项，不确定性集合不仅包含风电，还可以把碳价作为随机变量处理，此时模糊集的构造维度增加，但整体框架不需要改动。
• 与强化学习结合：把第二阶段问题替换为一个学习得到的决策策略网络，第一阶段仍然用优化求解，形成优化-学习的混合框架。这个方向在近期文献中很热，但其理论收敛性研究还不够成熟，建议先保证模型本身的数值稳定性再考虑引入学习模块。
• 更细粒度的热网动态建模：如果计算资源充足，可以尝试保留热网的二阶动态模型，但需要在求解前做好模型降阶或时空网格加密的自适应处理，否则求解时间会非常感人。

不要一次性把所有扩展都塞进一个模型里。先把基础的两阶段DRO框架彻底跑通，确认结果合理，再逐步增加复杂度，这样定位问题会容易得多。

我在实际调试中最深的一点体会是：分布鲁棒优化这类模型，真正难的往往不在理论推导，而在把模型转成可数值求解的实现细节——约束怎么写、变量怎么组织、求解器参数怎么调、对偶问题怎么验证，每一步都有可能让最终结果偏离预期。希望这篇内容能给你一条已经踩平的路，让你在Matlab代码实现时少走一些弯路。最后再分享一个小技巧：跑完主程序后，把每个场景下的二阶段调整成本单独导出来画个直方图，你会直观理解为什么单纯降低期望成本不一定能提升实际运行的稳定性。