近年来在做电热综合能源系统调度研究时我明显感觉到一个趋势纯物理建模的路子越来越走不动了。原因倒不复杂——热网动态特性、建筑热惯性、可再生出力波动这些环节实际运行数据里藏着大量模型难以精确刻画的信息尤其在样本量不足的小场景下传统数据驱动模型容易过拟合而纯物理模型又泛化不足。这篇要分享的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代码里我把第一阶段的决策变量定义为变量含义类型u_i机组i的启停状态0/1二进制x_i机组i的开机动作二进制H_CHPCHP机组的热出力计划连续H_TES_in / H_TES_out储热罐充/放热功率连续S_TES储热罐蓄热状态连续第一阶段的约束包括机组最小启停时间约束、CHP电热可行域约束、储热罐容量约束和热平衡约束。这些约束的共同特点是它们只依赖第一阶段的决策不涉及不确定参数的实现因此在优化开始前就可以确定。3.2 第二阶段决策变量与实时调整约束到了第二阶段风电出力ξ_s的真实值已经“揭晓”。此时系统需要在第一阶段决策的基础上对常规机组的电出力、CHP的电出力、电锅炉的用电功率等变量进行再调整。第二阶段的决策变量主要包括ΔP_i常规机组i的出力调整量ΔP_CHPCHP电出力调整量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是场景概率分布属于模糊集Ω。这个三层结构没法直接丢给求解器。标准处理步骤是把内层的max问题拆出来。由于内层是关于概率p的线性问题而模糊集由1-范数和∞-范数约束定义恰好是线性约束因此内层max问题可以等价变换为一个包含对偶变量的优化问题。把第二阶段的min Q(x, ξ_s)用其对偶问题替换与外层min合并。最终整个模型被转换为一个单层的混合整数线性规划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, T1, 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问题的对偶再将其并入主问题。我在代码里采用了列与约束生成CCG算法避免一次性把大规模问题直接塞给求解器。4.3 CCG主问题与子问题的迭代求解策略CCG算法的核心是“主问题-子问题”交替迭代主问题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 不同范数约束下的调度结果对比先看只使用单一范数约束的结果。模糊集类型第一阶段成本万元平均调整成本万元总成本万元相对确定性模型增幅确定性模型48.2048.2—1-范数θ0.551.37.658.922.2%∞-范数θ0.0552.88.961.728.0%1-范数∞-范数51.77.459.122.6%可以看到单一∞-范数的模型最保守总成本最高。组合约束的总成本只比单用1-范数高一点点但能有效限制单个场景的概率异常实际调度方案更稳健。从物理角度看这个结果也合理1-范数约束在控制整体分布偏移方面效率更高而∞-范数约束补上了它对单点概率失控的短板。两者组合能在大约相同总成本下提供更强的分布鲁棒性保证。5.3 迭代收敛曲线与求解耗时表现CCG算法在S100的场景数量下迭代9次收敛总耗时约186秒。每轮迭代的主问题规模虽然逐步增大但增加的都是线性约束和连续变量MILP部分增长有限。在实际运行中耗时最长的不是主问题而是子问题中第二阶段LP的批量求解。我对这部分做了场景并行化处理利用Matlab的parfor对场景循环并行计算子问题求解时间从原来的约15分钟压缩到不到3分钟。如果你也在做类似研究这一步优化非常值得投入时间。6. 实操中踩过的坑对偶转换与约束线性化的细节6.1 概率变量的非负约束丢失导致的对偶模型错误这个坑是我印象最深的。第一次做对偶转换时我手工推导后把概率变量的非负约束p ≥ 0给丢了导致生成的模糊集实际上允许负概率出现。症状表现为主问题迭代到某一步后子问题的目标值出现负的调整成本而且越迭代越离谱。排查半天才发现问题出在对偶转换时对偶变量符号方向搞反了。后来我长了个记性对偶变换后一定要先做小规模数值验证。具体做法是取一个S3的小算例手工列出原问题和转换后问题的KKT条件逐项检查约束对应关系是否完整。6.2 热网动态约束线性化过程中出现的数值病态电热综合能源系统的热网动态方程是偏微分方程在调度模型里通常采用节点法或有限差分法离散。问题在于当离散时间步长取的比较小时热网管道传输延迟会引入大量二进制变量导致模型求解时间爆炸。我采用的替代方案是把热网的传输延迟近似为多阶惯性环节的离散状态空间模型舍去热损非线性的高阶项只保留一次线性近似。这样处理之后模型规模显著下降而调度结果的偏差控制在2%以内在实际工程中可以接受。6.3 储热罐约束的时序耦合处理储热罐的状态变量在两个时间段之间存在耦合约束S_TES(t1) 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. 模型扩展思路从单区域向多区域与多能互补方向推进如果想把这套方法应用到你自己的课题中有几个扩展方向值得考虑多区域互联把单区域的电热耦合扩展为多个区域通过联络线和热网管道互联的情形模糊集的构造逻辑不变但需要额外处理区域间的耦合约束CCG子问题的规模会成倍增长。加入碳交易机制在目标函数中增加碳配额成本和碳交易成本项不确定性集合不仅包含风电还可以把碳价作为随机变量处理此时模糊集的构造维度增加但整体框架不需要改动。与强化学习结合把第二阶段问题替换为一个学习得到的决策策略网络第一阶段仍然用优化求解形成优化-学习的混合框架。这个方向在近期文献中很热但其理论收敛性研究还不够成熟建议先保证模型本身的数值稳定性再考虑引入学习模块。更细粒度的热网动态建模如果计算资源充足可以尝试保留热网的二阶动态模型但需要在求解前做好模型降阶或时空网格加密的自适应处理否则求解时间会非常感人。不要一次性把所有扩展都塞进一个模型里。先把基础的两阶段DRO框架彻底跑通确认结果合理再逐步增加复杂度这样定位问题会容易得多。我在实际调试中最深的一点体会是分布鲁棒优化这类模型真正难的往往不在理论推导而在把模型转成可数值求解的实现细节——约束怎么写、变量怎么组织、求解器参数怎么调、对偶问题怎么验证每一步都有可能让最终结果偏离预期。希望这篇内容能给你一条已经踩平的路让你在Matlab代码实现时少走一些弯路。最后再分享一个小技巧跑完主程序后把每个场景下的二阶段调整成本单独导出来画个直方图你会直观理解为什么单纯降低期望成本不一定能提升实际运行的稳定性。