干这个方向的人应该都有体会每次看到数据驱动两阶段分布鲁棒这几个词叠在一起第一反应都是这又能水一篇论文——但如果你真的动手用过Matlab去复现会发现里面有不少坑是论文里根本不写的。这篇博文我准备把这个课题从建模到求解再到代码实现完整拆开讲一遍围绕电热综合能源系统这个场景说明白1-范数和∞-范数约束在数据驱动分布鲁棒里面到底扮演什么角色以及你在Matlab里应当怎么落地。无论你是刚入门的研究生还是已经调过几版CCG算法的老手我相信这里都有你可以直接用的东西。数据驱动的两阶段分布鲁棒1-范数和∞-范数约束的电热综合能源系统研究Matlab代码实现1. 问题从哪来为什么电热综合能源系统需要分布鲁棒1.1 电热系统里的不确定性到底有多棘手先问一个问题一个电热综合能源系统Integrated Electricity and Heat System, IEHS里面最让调度员头疼的随机量是什么答案基本集中在两个——风电出力和热负荷。风电出力就不用细说了天然具有间歇性但热负荷的随机性容易被忽略它受天气、建筑保温、用户行为影响波动幅度甚至比电力负荷还大。而偏偏电和热在系统中是强耦合的联结二者的核心设备是热电联产机组CHP它的电出力和热出力在运行区间内互为约束不能单独调节。这意味着电源侧的风电不确定性、负荷侧的热需求波动会通过CHP的运行域传导到整个系统导致一步错、步步错。传统的处理手段有两类。一类是随机规划Stochastic Programming假设不确定量的概率分布是完全已知的另一类是鲁棒优化Robust Optimization只用一个不确定集合来描述所有可能取值完全放弃概率信息。前者的问题是真实系统中你根本拿不到精确分布顶多是从历史数据里估计一个经验分布。后者的问题是过于保守为了覆盖极端场景调度结果的经济性会变得很差而且这种保守是盲目的——它不区分这个场景概率很高和这个场景十年一遇。那有没有一种方法既利用数据里的概率信息又不完全依赖某个假设的分布这就是分布鲁棒优化Distributionally Robust Optimization, DRO登场的理由。它假设不确定量的真实分布属于一个模糊集Ambiguity Set这个集合由历史数据构造覆盖所有可能的分布。我们优化的目标是在最坏情况分布下让总成本最小。换句话说DRO站在随机规划和鲁棒优化中间——用数据构造模糊集在模糊集内寻找最坏分布从而兼顾鲁棒性和经济性。1.2 为什么偏偏选两阶段再来看两阶段三个字。在电热系统的日前调度中决策天然分两个阶段第一阶段日前决策在不确定性实现之前决定机组启停状态、各机组的基础出力基准点、储热罐的充放热计划。这些决策一旦做出当天之内基本无法调整。第二阶段实时调整风电出力、热负荷的实际值逐渐明确后系统通过调整机组出力偏差、切负荷、弃风等手段保证供需平衡。这个阶段的成本是情景依赖的不同场景下的调整代价不同。用一个直白类比第一阶段的决策好比你提前订好了餐厅和菜品不能退了第二阶段是你落座后根据当天胃口好坏决定加菜还是打包可以灵活应对。好的调度策略就是让第一阶段决策尽量平衡使得第二阶段无论出现什么情况总代价都可控。两阶段分布鲁棒模型就是在第一阶段就考虑第二阶段的最坏情况期望成本从而导出在不确定性下长期表现最优的日前计划。这两个阶段合在一起正好对应了一个整数变量加连续变量的混合整数线性规划MILP问题——而这也是为什么这个项目能在Matlab中用YALMIP配合求解器来解决的原因。2. 数学模型拆解1-范数和∞-范数约束的真实含义2.1 模糊集怎么构造才叫数据驱动正式写模型之前必须把一个概念搞清楚分布鲁棒优化里的数据驱动本质上是说模糊集不是拍脑袋定的而是从历史数据中生成的。第一步把风电出力和热负荷的历史数据整理成N个离散场景每个场景对应一组不确定参数取值。假设每个场景的权重即概率值为(p_k)初始经验分布为(p_0)其中(p_{0,k} 1/N)。然后关键来了真实分布不一定等于经验分布但我们相信它不会离经验分布太远。怎么度量太远论文里常用的是概率距离或范数距离。这个课题选用的是1-范数和∞-范数同时约束模糊集可以写成[ D \left{ p \in \mathbb{R}^{N} \ \middle| \ | p - p_0 |1 \le \theta_1, \ | p - p_0 |{\infty} \le \theta_{\infty}, \ \sum_{k1}^{N} p_k 1, \ p_k \ge 0 \right} ]其中(| p - p_0 |1 \sum{k1}^{N} |p_k - p_{0,k}|)约束所有场景概率偏差的绝对值之和(| p - p_0 |{\infty} \max{k} |p_k - p_{0,k}|)约束单个场景概率偏差的最大值(\theta_1) 和 (\theta_{\infty}) 分别是两个范数约束的半径允许偏差上界。这个构造方式不是随便选的它有两个非常实际的好处。第一相比Wasserstein距离模糊集或KL散度模糊集范数约束模糊集的可行域是一个多面体线性规划求解方便第二1-范数和∞-范数是互补的——1-范数约束限制整体偏移总量∞-范数约束防止某个点的概率被单独拉得太极端。两者结合模糊集既不宽得离谱又不会出现某个场景概率被调成0或者过高的病态情况。2.2 两阶段模型的完整表达式有了模糊集两阶段分布鲁棒模型可以写成如下形式[ \min_{x \in X} \ c^T x \max_{p \in D} \sum_{k1}^{N} p_k \cdot Q(x, \xi_k) ]其中(x) 表示第一阶段决策变量(X) 是满足机组启停、爬坡、储热罐容量等约束的可行域(c^T x) 是第一阶段的燃料成本、启停成本和购电成本(Q(x, \xi_k)) 是在第一阶段决策固定为x、不确定量取第k个场景值(\xi_k)时第二阶段的最优运行调整成本。这个式子的核心逻辑是外层先选第一阶段计划内层在最坏分布下评估期望调整成本。特别要注意(\max_{p \in D})把分布概率作为决策变量因此这里的期望不是对单一固定分布取期望而是在一系列可能分布中选那个让成本最大的分布来计算期望——这正是DRO保守性的来源也是它相比随机规划能给出稳妥方案的原因。再展开第二阶段问题(Q(x, \xi_k))它是一个线性规划目标函数为[ Q(x, \xi_k) \min_y \ d^T y ]约束条件包括电力平衡约束(\sum_{\text{unit}} P_{\text{unit}} P_{\text{wind},k} - P_{\text{curtail},k} P_{\text{load},k} P_{\text{HP}} P_{\text{EB}})热力平衡约束(\sum_{\text{CHP}} H_{\text{CHP}} H_{\text{GB}} H_{\text{discharge},k} - H_{\text{charge},k} H_{\text{load},k})机组出力调整上下限、CHP电热运行可行域、储能约束等其中(y)是第二阶段的调整决策变量各机组出力偏差、弃风量、切负荷量、储热罐充放热功率(d)是相应的惩罚系数比如弃风惩罚、切负荷惩罚通常设置得远高于正常发电成本以保证调度机制不乱切。这里还有值得一提的细节第一阶段和第二阶段通过机组出力基准点和储热罐充放热计划耦合。如果第一阶段定的储热罐放热计划不合理第二阶段热负荷一波动要么切负荷要么启停燃气锅炉成本瞬间飙升。DRO的价值恰恰在这里——它会逼着第一阶段在模糊集内把所有场景都想一遍得到一个无论哪个场景真实发生都不至于太狼狈的储热计划和出力计划。2.3 为什么是1-范数和∞-范数而非其它很多人第一次看文献会问为什么不用Wasserstein距离为什么不用矩模糊集选择1-范数和∞-范数有几个实际原因线性性这两个范数约束在多面体上定义线性规划或MILP求解器可以直接处理不需要引入半定规划或二次约束。对于系统规模较大的电热综合能源系统这关系到能不能在可接受时间内求得解。对偶友好在CCG算法框架下内层求最坏分布是一个线性规划其对偶问题的结构非常清晰可以用对偶变量把max-min问题化成单层可解形式。直观可解释(\theta_1)直接告诉你所有场景真实概率相对经验概率的总偏差上界(\theta_\infty)告诉你在单个场景上最大偏差的容忍度。这比纯抽象的距离度量更容易面向工程人员解释和调参。当然它的代价是模糊集的形状比较方不如Wasserstein球那样在连续概率空间里有漂亮的度量性质。但如果你做的是电热系统日前调度这类工程落地型课题实用性和可求解性往往比纯粹的理论优雅更重要——这也是我在复现这个方向时觉得范数约束非常适合作为入门又足够支撑一篇论文的原因。3. 求解算法核心CCG迭代与Matlab实现细节3.1 先把max-min问题重塑成可以解的形态模型写出来是一回事能求解是另一回事。直接面对上面的两阶段模型内层(\max_{p \in D} \sum p_k Q(x, \xi_k))本质上是个嵌套优化问题不能直接丢给Gurobi或CPLEX解。标准解法是列与约束生成Column-and-Constraint Generation, CCG它的思路和Benders分解有亲缘关系但收敛速度通常更快。CCG的核心思想是把原问题拆成一个主问题MP和一个子问题SP。主问题包含第一阶段决策变量和一个有限的最坏场景候选集合每次迭代后把SP找出的最坏场景作为一个新的约束或列加入主问题子问题给定第一阶段解寻找使第二阶段成本最大的场景概率分布并把该分布情形反馈给主问题。迭代收敛判据是主问题的目标函数值下界与子问题得到的真实第二阶段最坏期望成本之和上界之差小于给定间隙(\epsilon)。具体到本项目子问题给定x后为[ \max_{p \in D} \sum_{k1}^N p_k \cdot Q(x, \xi_k) ]其中(Q(x, \xi_k))本身是线性规划的最优值函数。通过线性规划强对偶可以把每个(Q(x, \xi_k))改写为最大化或最小化的对偶形式然后与外层的(\max_{p \in D})合并成一个更大的线性规划。这样SP就变成一个标准LP可以直接用求解器解出最坏分布(p^)。随后MP在下一轮迭代中加入与(p^)对应的约束[ c^T x \sum_{k1}^{N} p_k^* \cdot Q(x, \xi_k) \le \eta ]并继续求解直到满足收敛判据。这个主问题加一列约束、子问题解一个LP的结构在Matlab中用YALMIP建模非常自然。3.2 CCG算法的逐步流程与Matlab代码骨架下面给出一个贴近实际工程的CCG迭代流程并配上Matlab核心代码骨架。Step 1初始化生成历史数据场景集({\xi_k}_{k1}^N)设置模糊集半径(\theta_1, \theta_\infty)初始化上界(UB\infty)、下界(LB-\infty)、迭代次数(iter1)设置收敛间隙(\epsilon)比如1e-4。Step 2求解主问题MP主问题是在当前最坏场景集合下求最低成本并把辅助变量(\eta)加入目标函数[ \min_{x, y, \eta} \ c^T x \eta ]主问题的YALMIP建模如下% 主问题变量 x sdpvar(ng, 1); % 第一阶段机组基准出力 u binvar(ng, 1); % 机组启停状态 eta sdpvar(1, 1); % 第二阶段成本辅助变量 % 目标函数 obj_MP c * x sum(StartCost .* u) eta; % 基本约束 Constraints_MP [A_eq * x b_eq, A_ineq * x b_ineq]; % CCG迭代中添加的最坏场景约束 % 每轮CCG迭代会追加一组约束下详 ... ops sdpsettings(solver, cplex, verbose, 1, display, iter); optimize(Constraints_MP, obj_MP, ops);Step 3求解子问题SP子问题解出最坏概率分布p。这里为了求解方便把第二阶段问题的对偶式加入将max-min转成LP% 子问题给定第一阶段解x_hat求最坏分布 p sdpvar(N, 1); % 场景概率变量 lambda sdpvar(N, 1); % 对偶变量每个场景对应一个 % 第二阶段成本Q_k用对偶表示Q_k max over dual variables % 这里示意性地将Q_k表示成对偶目标值的线性形式 Q_k dual_expr; % 由具体模型推导 % 内层 worst-case 期望 obj_SP - sum(p .* Q_k); % 若Q_k是对偶最大化则这里是凸的 % 模糊集约束1-范数 ∞-范数 Constraints_SP [sum(p) 1, p 0, ...]; Constraints_SP [Constraints_SP, norm(p - p0, 1) theta1]; Constraints_SP [Constraints_SP, norm(p - p0, inf) theta_inf]; optimize(Constraints_SP, obj_SP, ops); p_star value(p);Step 4更新上界用当前第一阶段解(x^*)直接计算所有场景下第二阶段真实成本结合第一阶段成本构成上界[ UB \min(UB, \ c^T x^* \sum_{k1}^N p_k^* \cdot Q(x^*, \xi_k)) ]Step 5更新下界与收敛判断主问题最优值就是当前下界(LB \text{obj_MP})。如果(UB - LB \le \epsilon)则停止否则把子问题找到的最坏分布对应的约束加入主问题就是给(p^*)固定后把第二阶段成本表达式作为约束加入MP迭代次数加1回到Step 2。3.3 代码里最容易踩的坑这个算法流程看起来不算复杂但真正在Matlab里把代码调通有几个坑我踩过很多次值得专门提醒。坑一第二阶段问题不可行。求解子问题时会发现某些场景下第一阶段给定的调度计划根本满足不了平衡约束。这通常是因为第一阶段忘记预留足够的调节空间或者机组出力区间设置得不合理。解决办法是在第二阶段目标中加入松弛变量切负荷量、弃风量并设一个很大的惩罚系数。这不仅是数值技巧在工程上也有实际意义——紧急情况下切一点负荷总比全系统崩溃要好。坑二YALMIP的norm函数在大规模场景下会生成很多额外变量。1-范数和∞-范数可以用线性约束等价替换不一定要用norm函数。比如1-范数约束等价于引入辅助变量(t_k\ge0)加上(p_k - p_{0,k} \le t_k)和(p_{0,k} - p_k \le t_k)以及(\sum t_k \le \theta_1)。用显式线性不等式实现求解速度有时比直接调norm快不少尤其当N达到几百上千的时候。坑三CCG初始迭代阶段上下界振荡迟迟不收敛。常见的原因是子问题求解精度不够或者主问题加入的约束写成了严格不等式。另一个原因是模糊集半径设置偏大导致子问题反复找极端场景主问题不断加约束迭代步数暴涨。实践中我会先检查半径是否合理再检查约束是否添加正确。坑四YALMIP CPLEX/Gurobi的版本兼容问题。尤其是2023版之后一些老代码里调用的sdpsettings参数名可能失效。建议在代码开头统一设置求解器为cplex或gurobi并把输出信息重定向到文件方便回溯调试。4. 系统建模电热综合能源系统的约束与参数怎么定4.1 系统组成与设备运行域算例系统不宜太复杂适合复现和迭代调试。我通常采用改进的IEEE 6节点电力系统耦合一个4节点热力系统。主要设备包括两台常规燃煤机组G1、G2作为第一阶段主力调节电源一台抽凝式CHP机组电热联合运行可行域是典型的梯形区域一台燃气锅炉GB作为热力备用一个储热罐HS承担热负荷转移功能风电场WF输出功率为不确定量电力负荷和热负荷其中热负荷也带有一定随机性。关键设备运行约束要写对尤其是CHP的运行域。抽凝式CHP的可行电出力与热出力关系可以近似表示为[ P_{\text{CHP}} \in [\underline{P}{\text{CHP}}(H{\text{CHP}}),\ \overline{P}{\text{CHP}}(H{\text{CHP}})] ]即电出力上下界都是热出力的线性函数构成一个梯形。这组约束决定了系统的电热耦合强度——CHP在满足热负荷的同时实际上绑架了部分电出力这是电热综合能源系统调度难度的根源。4.2 不确定性数据生成与场景缩减要体现数据驱动首先要有历史数据场景集。实际项目里我们可能只有预测曲线和误差统计信息因此常通过蒙特卡洛或拉丁超立方采样生成N个等概率样本。例如风电出力可表示为预测值加误差[ P_{\text{wind},k} P_{\text{wind,forecast}} \varepsilon_k,\quad \varepsilon_k \sim \mathcal{N}(0, \sigma_w^2) ]热负荷同样处理。要注意样本数量N不能太小——范数约束模糊集的经验分布是以N个场景为支撑的N太小则经验分布可信度低模糊集半径也难以合理化。我一般取N200~500。如果N太大第二阶段LP的求解规模会上升不过现在CPLEX解几百个场景的LP是没问题的影响不大。场景生成好后对每个场景做归一化确保数值稳定性。这是个大坑——不归一化的话热负荷数值可能上千风电出力数值却只有几十两者一起进入模糊集的范数距离计算时风电的偏差会被热负荷的绝对值掩盖导致模糊集形同虚设。4.3 存储与参数标定给出一个实用的参数配置表可以直接用于复现。这个表是我在调参中积累的默认档位适合做一个基准算例后续再根据你的研究需要调整。参数名称数值说明场景数 N200历史数据采样数θ1 (1-范数半径)0.6总概率偏差上限θ∞ (∞-范数半径)0.05单场景概率偏差上限收敛间隙 ε1e-4CCG迭代停止条件弃风惩罚系数50 $/MWh第二阶段目标惩罚切负荷惩罚系数500 $/MWh远高于正常成本防止随意切负荷储热罐容量200 MWh储热罐可存储热量上限CHP最大电出力200 MW—燃气锅炉最大热出力150 MWth—注意θ1和θ∞不是拍脑袋定的需要通过置信水平反推。常用的经验公式是在置信水平α下θ1满足(P{|p-p_0|_1 \le \theta_1} \ge \alpha)。这需要用到概率不等式估算。工程上最简单的方式是换几组θ值跑一遍观察结果成本的灵敏度。θ越大越保守成本越高选在成本急剧增加之前的拐点附近通常是一个合理选择。这一点我会在后面详细展开。5. 算例结果解读怎么证明分布鲁棒比随机规划更好5.1 对比实验设计想让论文或报告有说服力就要做对比。通常设置三组方案随机规划SP直接用经验分布(p_0)作为唯一分布不走模糊集传统鲁棒优化RO不确定量落在区间内不区分概率最坏情况下取极端场景两阶段分布鲁棒DRO本文所用模型模糊集由1-范数和∞-范数约束构造。每组方案跑完后对比三个指标总成本、第一阶段计划、第二阶段最坏期望成本。实际算例跑出来的典型结果是SP的总成本最低但它的风险也最大——一旦真实分布偏离经验分布实际运行成本会急剧上升RO最保守总成本明显偏高第一阶段几乎为所有可能场景都预留了大量备用导致正常日子浪费DRO介于两者之间成本只比SP高几个百分点却能保证在模糊集范围内不出现成本失控。这个结果其实是符合直觉的DRO用数据的概率信息换取了鲁棒性而且通过调节θ1和θ∞可以连续地在SP和RO之间移动——θ越小方案越接近SPθ越大方案越接近RO。这种可调性在工程决策中非常宝贵。5.2 置信度分析为什么不真实分布怎么应对再做一个更有说服力的验证通过蒙特卡洛模拟随机生成大量真实分布不一定等于经验分布然后分别把SP、RO、DRO三套方案放到这些真实分布下计算实际运行成本统计成本分布。结果一般会显示SP方案的成本分布尾部很长——偶尔出现极端高成本RO方案的成本分布集中但均值很高DRO方案的成本分布尾部截断明显优于SP均值又显著低于RO。这种尾部风险控制正是分布鲁棒真正的卖点不是压最低平均成本而是在不严重牺牲经济性的前提下控制最坏情况下的实际损失。如果你的算例结果呈现这个趋势说明你的模型和代码写对了。5.3 参数灵敏度分析θ怎么选才合理θ1和θ∞的选取没有绝对标准但做好灵敏度分析是论文必备内容。我也是在反复实验中发现θ1对结果的影响通常比θ∞更明显因为1-范数约束控制整体概率偏移的总量而∞-范数主要限制单个场景的极端偏移。实际操作时可以固定θ∞不变把θ1从0.2、0.4、0.6、0.8、1.0依次增大观察总成本和迭代次数的变化。典型趋势是θ1增大到某个阈值之前总成本缓慢上升超过某个阈值后成本可能大幅跳增此时模糊集过大、方案过度保守。工程上选择阈值之前的拐点作为推荐值比如0.6。同时θ1增大会让CCG迭代次数增加因为子问题更容易找出不同的极端场景。这个成本-鲁棒性-计算量的综合权衡需要在报告里用图表呈现结论才扎实。6. Matlab代码调试中的常见问题与排查技巧6.1 问题速查表现象可能原因解决办法子问题LP求解失败或不可行模糊集为空θ1太小或θ∞太小导致p无法满足所有约束增大θ1或θ∞检查sum(p)1与范数约束是否冲突CCG迭代次数过多50模糊集半径过大场景数量过多收敛间隙太小适当减小θ1对场景数做缩减放宽ε到1e-3上界持续不下降主问题没有加入正确的割约束第二阶段成本计算与子问题不一致核对主问题新增约束是否使用了子问题返回的p_star检查Q_k表达式是否一致求解结果出现负成本或异常大数目标函数符号写反惩罚系数设置过大导致数值问题逐项检查目标函数调整惩罚系数到合理范围YALMIP报Unable to instantiate solver求解器路径未配置或版本不兼容在yalmip/operators中检查求解器是否支持LP/MILP重新添加CPLEX/Gurobi路径6.2 我调试了几个小时才发现的细节有几个非常不起眼的细节浪费过我大量时间这里专门写出来。第一个是关于YALMIP中repmat的使用。当第二阶段约束里涉及多个场景的表达式时很多人习惯用循环逐个写约束。场景数少还行场景数一多模型构建时间成倍增长。更好的做法是向量化用repmat或直接构造二维变量矩阵一次性写入所有场景的约束。这么做模型构建时间可以从几十秒降到一两秒。第二个是通过value()函数断言检查结果。每次迭代后手动校验一下场景概率p_star是否满足sum(p)1是否都在[0,1]区间内。这个习惯帮我发现过两次模糊集约束写错的问题——一次是因为缺少sum(p)1另一次是因为∞-范数约束误写成了逐元素约束的平方和。第三个是关于求解器日志。CPLEX在求解MILP时日志非常长但关键信息是Best bound和Best integer两个值。建议写一个简单的日志解析脚本自动提取这两个数并绘制收敛曲线能直观判断CCG和MILP的收敛趋势。没有这条曲线你很难判断迭代是在收敛还是卡住了。第四个是一个容易混淆的点1-范数约束中的θ1与∞-范数中的θ∞它们与场景数的关系需要单独校准。当N改变时经验分布p0的各分量p0_k1/N也会变所以同样的θ1和θ∞在不同N下对应不同的模糊程度。做灵敏度分析时务必固定N否则结果不可比。6.3 提升代码运行效率的三个小技巧场景缩减如果初始采样场景数达到1000可以先做k-means聚类或同步回代缩减SBR把场景数降到200~300再进入模糊集模型。这样CCG的迭代次数和子问题规模都会显著降低。热启动Warm StartCCG每次迭代求解主问题时可以以上一轮的解作为初始点YALMIP的assign和optimize的usex0参数可以配合使用。在MILP问题中热启动效果尤其明显。缓存子问题的对偶解在CCG迭代过程中相同场景在相邻迭代中可能被反复求解如果能够缓存每个场景的对偶变量初始值可以加速LP求解。虽然CPLEX自带了一些记忆机制但在自定义实现里显式复用初始点仍然能带来可感知的提升。7. 我个人的几点实操体会这个项目做下来我最大的体会是数据驱动两阶段分布鲁棒优化并不像有些论文写得那么玄乎它本质上就是用数据约束概率分布的不确定性范围然后在这个范围内做最坏情况决策。它的难点反而在细节——模糊集怎么表达、对偶怎么转、CCG的上下界怎么算、YALMIP的向量化怎么写每一步都有不止一个坑。但把这些坑踩平之后你会发现这套框架在电热综合能源系统里特别好用因为CHP的强耦合特性和储热装置的时间转移特性恰恰能让两阶段结构发挥出最大价值。最后再分享一个小技巧当你面对自己的算例结果不知道如何解释时先不要急着调参数而是用蒙特卡洛法把SP、RO、DRO三套方案放到真实分布下做一次批量模拟画一张成本分布箱线图。这张图基本能一眼看出DRO的价值也是论文中最容易打动审稿人的一张图。很多时候方案好坏的差异并不在第一阶段的计划表上而在于它面对真实不确定性时的耐受力——这个道理做调度的人应该都懂。
数据驱动两阶段分布鲁棒优化在电热综合能源系统中的应用及Matlab实现
干这个方向的人应该都有体会每次看到数据驱动两阶段分布鲁棒这几个词叠在一起第一反应都是这又能水一篇论文——但如果你真的动手用过Matlab去复现会发现里面有不少坑是论文里根本不写的。这篇博文我准备把这个课题从建模到求解再到代码实现完整拆开讲一遍围绕电热综合能源系统这个场景说明白1-范数和∞-范数约束在数据驱动分布鲁棒里面到底扮演什么角色以及你在Matlab里应当怎么落地。无论你是刚入门的研究生还是已经调过几版CCG算法的老手我相信这里都有你可以直接用的东西。数据驱动的两阶段分布鲁棒1-范数和∞-范数约束的电热综合能源系统研究Matlab代码实现1. 问题从哪来为什么电热综合能源系统需要分布鲁棒1.1 电热系统里的不确定性到底有多棘手先问一个问题一个电热综合能源系统Integrated Electricity and Heat System, IEHS里面最让调度员头疼的随机量是什么答案基本集中在两个——风电出力和热负荷。风电出力就不用细说了天然具有间歇性但热负荷的随机性容易被忽略它受天气、建筑保温、用户行为影响波动幅度甚至比电力负荷还大。而偏偏电和热在系统中是强耦合的联结二者的核心设备是热电联产机组CHP它的电出力和热出力在运行区间内互为约束不能单独调节。这意味着电源侧的风电不确定性、负荷侧的热需求波动会通过CHP的运行域传导到整个系统导致一步错、步步错。传统的处理手段有两类。一类是随机规划Stochastic Programming假设不确定量的概率分布是完全已知的另一类是鲁棒优化Robust Optimization只用一个不确定集合来描述所有可能取值完全放弃概率信息。前者的问题是真实系统中你根本拿不到精确分布顶多是从历史数据里估计一个经验分布。后者的问题是过于保守为了覆盖极端场景调度结果的经济性会变得很差而且这种保守是盲目的——它不区分这个场景概率很高和这个场景十年一遇。那有没有一种方法既利用数据里的概率信息又不完全依赖某个假设的分布这就是分布鲁棒优化Distributionally Robust Optimization, DRO登场的理由。它假设不确定量的真实分布属于一个模糊集Ambiguity Set这个集合由历史数据构造覆盖所有可能的分布。我们优化的目标是在最坏情况分布下让总成本最小。换句话说DRO站在随机规划和鲁棒优化中间——用数据构造模糊集在模糊集内寻找最坏分布从而兼顾鲁棒性和经济性。1.2 为什么偏偏选两阶段再来看两阶段三个字。在电热系统的日前调度中决策天然分两个阶段第一阶段日前决策在不确定性实现之前决定机组启停状态、各机组的基础出力基准点、储热罐的充放热计划。这些决策一旦做出当天之内基本无法调整。第二阶段实时调整风电出力、热负荷的实际值逐渐明确后系统通过调整机组出力偏差、切负荷、弃风等手段保证供需平衡。这个阶段的成本是情景依赖的不同场景下的调整代价不同。用一个直白类比第一阶段的决策好比你提前订好了餐厅和菜品不能退了第二阶段是你落座后根据当天胃口好坏决定加菜还是打包可以灵活应对。好的调度策略就是让第一阶段决策尽量平衡使得第二阶段无论出现什么情况总代价都可控。两阶段分布鲁棒模型就是在第一阶段就考虑第二阶段的最坏情况期望成本从而导出在不确定性下长期表现最优的日前计划。这两个阶段合在一起正好对应了一个整数变量加连续变量的混合整数线性规划MILP问题——而这也是为什么这个项目能在Matlab中用YALMIP配合求解器来解决的原因。2. 数学模型拆解1-范数和∞-范数约束的真实含义2.1 模糊集怎么构造才叫数据驱动正式写模型之前必须把一个概念搞清楚分布鲁棒优化里的数据驱动本质上是说模糊集不是拍脑袋定的而是从历史数据中生成的。第一步把风电出力和热负荷的历史数据整理成N个离散场景每个场景对应一组不确定参数取值。假设每个场景的权重即概率值为(p_k)初始经验分布为(p_0)其中(p_{0,k} 1/N)。然后关键来了真实分布不一定等于经验分布但我们相信它不会离经验分布太远。怎么度量太远论文里常用的是概率距离或范数距离。这个课题选用的是1-范数和∞-范数同时约束模糊集可以写成[ D \left{ p \in \mathbb{R}^{N} \ \middle| \ | p - p_0 |1 \le \theta_1, \ | p - p_0 |{\infty} \le \theta_{\infty}, \ \sum_{k1}^{N} p_k 1, \ p_k \ge 0 \right} ]其中(| p - p_0 |1 \sum{k1}^{N} |p_k - p_{0,k}|)约束所有场景概率偏差的绝对值之和(| p - p_0 |{\infty} \max{k} |p_k - p_{0,k}|)约束单个场景概率偏差的最大值(\theta_1) 和 (\theta_{\infty}) 分别是两个范数约束的半径允许偏差上界。这个构造方式不是随便选的它有两个非常实际的好处。第一相比Wasserstein距离模糊集或KL散度模糊集范数约束模糊集的可行域是一个多面体线性规划求解方便第二1-范数和∞-范数是互补的——1-范数约束限制整体偏移总量∞-范数约束防止某个点的概率被单独拉得太极端。两者结合模糊集既不宽得离谱又不会出现某个场景概率被调成0或者过高的病态情况。2.2 两阶段模型的完整表达式有了模糊集两阶段分布鲁棒模型可以写成如下形式[ \min_{x \in X} \ c^T x \max_{p \in D} \sum_{k1}^{N} p_k \cdot Q(x, \xi_k) ]其中(x) 表示第一阶段决策变量(X) 是满足机组启停、爬坡、储热罐容量等约束的可行域(c^T x) 是第一阶段的燃料成本、启停成本和购电成本(Q(x, \xi_k)) 是在第一阶段决策固定为x、不确定量取第k个场景值(\xi_k)时第二阶段的最优运行调整成本。这个式子的核心逻辑是外层先选第一阶段计划内层在最坏分布下评估期望调整成本。特别要注意(\max_{p \in D})把分布概率作为决策变量因此这里的期望不是对单一固定分布取期望而是在一系列可能分布中选那个让成本最大的分布来计算期望——这正是DRO保守性的来源也是它相比随机规划能给出稳妥方案的原因。再展开第二阶段问题(Q(x, \xi_k))它是一个线性规划目标函数为[ Q(x, \xi_k) \min_y \ d^T y ]约束条件包括电力平衡约束(\sum_{\text{unit}} P_{\text{unit}} P_{\text{wind},k} - P_{\text{curtail},k} P_{\text{load},k} P_{\text{HP}} P_{\text{EB}})热力平衡约束(\sum_{\text{CHP}} H_{\text{CHP}} H_{\text{GB}} H_{\text{discharge},k} - H_{\text{charge},k} H_{\text{load},k})机组出力调整上下限、CHP电热运行可行域、储能约束等其中(y)是第二阶段的调整决策变量各机组出力偏差、弃风量、切负荷量、储热罐充放热功率(d)是相应的惩罚系数比如弃风惩罚、切负荷惩罚通常设置得远高于正常发电成本以保证调度机制不乱切。这里还有值得一提的细节第一阶段和第二阶段通过机组出力基准点和储热罐充放热计划耦合。如果第一阶段定的储热罐放热计划不合理第二阶段热负荷一波动要么切负荷要么启停燃气锅炉成本瞬间飙升。DRO的价值恰恰在这里——它会逼着第一阶段在模糊集内把所有场景都想一遍得到一个无论哪个场景真实发生都不至于太狼狈的储热计划和出力计划。2.3 为什么是1-范数和∞-范数而非其它很多人第一次看文献会问为什么不用Wasserstein距离为什么不用矩模糊集选择1-范数和∞-范数有几个实际原因线性性这两个范数约束在多面体上定义线性规划或MILP求解器可以直接处理不需要引入半定规划或二次约束。对于系统规模较大的电热综合能源系统这关系到能不能在可接受时间内求得解。对偶友好在CCG算法框架下内层求最坏分布是一个线性规划其对偶问题的结构非常清晰可以用对偶变量把max-min问题化成单层可解形式。直观可解释(\theta_1)直接告诉你所有场景真实概率相对经验概率的总偏差上界(\theta_\infty)告诉你在单个场景上最大偏差的容忍度。这比纯抽象的距离度量更容易面向工程人员解释和调参。当然它的代价是模糊集的形状比较方不如Wasserstein球那样在连续概率空间里有漂亮的度量性质。但如果你做的是电热系统日前调度这类工程落地型课题实用性和可求解性往往比纯粹的理论优雅更重要——这也是我在复现这个方向时觉得范数约束非常适合作为入门又足够支撑一篇论文的原因。3. 求解算法核心CCG迭代与Matlab实现细节3.1 先把max-min问题重塑成可以解的形态模型写出来是一回事能求解是另一回事。直接面对上面的两阶段模型内层(\max_{p \in D} \sum p_k Q(x, \xi_k))本质上是个嵌套优化问题不能直接丢给Gurobi或CPLEX解。标准解法是列与约束生成Column-and-Constraint Generation, CCG它的思路和Benders分解有亲缘关系但收敛速度通常更快。CCG的核心思想是把原问题拆成一个主问题MP和一个子问题SP。主问题包含第一阶段决策变量和一个有限的最坏场景候选集合每次迭代后把SP找出的最坏场景作为一个新的约束或列加入主问题子问题给定第一阶段解寻找使第二阶段成本最大的场景概率分布并把该分布情形反馈给主问题。迭代收敛判据是主问题的目标函数值下界与子问题得到的真实第二阶段最坏期望成本之和上界之差小于给定间隙(\epsilon)。具体到本项目子问题给定x后为[ \max_{p \in D} \sum_{k1}^N p_k \cdot Q(x, \xi_k) ]其中(Q(x, \xi_k))本身是线性规划的最优值函数。通过线性规划强对偶可以把每个(Q(x, \xi_k))改写为最大化或最小化的对偶形式然后与外层的(\max_{p \in D})合并成一个更大的线性规划。这样SP就变成一个标准LP可以直接用求解器解出最坏分布(p^)。随后MP在下一轮迭代中加入与(p^)对应的约束[ c^T x \sum_{k1}^{N} p_k^* \cdot Q(x, \xi_k) \le \eta ]并继续求解直到满足收敛判据。这个主问题加一列约束、子问题解一个LP的结构在Matlab中用YALMIP建模非常自然。3.2 CCG算法的逐步流程与Matlab代码骨架下面给出一个贴近实际工程的CCG迭代流程并配上Matlab核心代码骨架。Step 1初始化生成历史数据场景集({\xi_k}_{k1}^N)设置模糊集半径(\theta_1, \theta_\infty)初始化上界(UB\infty)、下界(LB-\infty)、迭代次数(iter1)设置收敛间隙(\epsilon)比如1e-4。Step 2求解主问题MP主问题是在当前最坏场景集合下求最低成本并把辅助变量(\eta)加入目标函数[ \min_{x, y, \eta} \ c^T x \eta ]主问题的YALMIP建模如下% 主问题变量 x sdpvar(ng, 1); % 第一阶段机组基准出力 u binvar(ng, 1); % 机组启停状态 eta sdpvar(1, 1); % 第二阶段成本辅助变量 % 目标函数 obj_MP c * x sum(StartCost .* u) eta; % 基本约束 Constraints_MP [A_eq * x b_eq, A_ineq * x b_ineq]; % CCG迭代中添加的最坏场景约束 % 每轮CCG迭代会追加一组约束下详 ... ops sdpsettings(solver, cplex, verbose, 1, display, iter); optimize(Constraints_MP, obj_MP, ops);Step 3求解子问题SP子问题解出最坏概率分布p。这里为了求解方便把第二阶段问题的对偶式加入将max-min转成LP% 子问题给定第一阶段解x_hat求最坏分布 p sdpvar(N, 1); % 场景概率变量 lambda sdpvar(N, 1); % 对偶变量每个场景对应一个 % 第二阶段成本Q_k用对偶表示Q_k max over dual variables % 这里示意性地将Q_k表示成对偶目标值的线性形式 Q_k dual_expr; % 由具体模型推导 % 内层 worst-case 期望 obj_SP - sum(p .* Q_k); % 若Q_k是对偶最大化则这里是凸的 % 模糊集约束1-范数 ∞-范数 Constraints_SP [sum(p) 1, p 0, ...]; Constraints_SP [Constraints_SP, norm(p - p0, 1) theta1]; Constraints_SP [Constraints_SP, norm(p - p0, inf) theta_inf]; optimize(Constraints_SP, obj_SP, ops); p_star value(p);Step 4更新上界用当前第一阶段解(x^*)直接计算所有场景下第二阶段真实成本结合第一阶段成本构成上界[ UB \min(UB, \ c^T x^* \sum_{k1}^N p_k^* \cdot Q(x^*, \xi_k)) ]Step 5更新下界与收敛判断主问题最优值就是当前下界(LB \text{obj_MP})。如果(UB - LB \le \epsilon)则停止否则把子问题找到的最坏分布对应的约束加入主问题就是给(p^*)固定后把第二阶段成本表达式作为约束加入MP迭代次数加1回到Step 2。3.3 代码里最容易踩的坑这个算法流程看起来不算复杂但真正在Matlab里把代码调通有几个坑我踩过很多次值得专门提醒。坑一第二阶段问题不可行。求解子问题时会发现某些场景下第一阶段给定的调度计划根本满足不了平衡约束。这通常是因为第一阶段忘记预留足够的调节空间或者机组出力区间设置得不合理。解决办法是在第二阶段目标中加入松弛变量切负荷量、弃风量并设一个很大的惩罚系数。这不仅是数值技巧在工程上也有实际意义——紧急情况下切一点负荷总比全系统崩溃要好。坑二YALMIP的norm函数在大规模场景下会生成很多额外变量。1-范数和∞-范数可以用线性约束等价替换不一定要用norm函数。比如1-范数约束等价于引入辅助变量(t_k\ge0)加上(p_k - p_{0,k} \le t_k)和(p_{0,k} - p_k \le t_k)以及(\sum t_k \le \theta_1)。用显式线性不等式实现求解速度有时比直接调norm快不少尤其当N达到几百上千的时候。坑三CCG初始迭代阶段上下界振荡迟迟不收敛。常见的原因是子问题求解精度不够或者主问题加入的约束写成了严格不等式。另一个原因是模糊集半径设置偏大导致子问题反复找极端场景主问题不断加约束迭代步数暴涨。实践中我会先检查半径是否合理再检查约束是否添加正确。坑四YALMIP CPLEX/Gurobi的版本兼容问题。尤其是2023版之后一些老代码里调用的sdpsettings参数名可能失效。建议在代码开头统一设置求解器为cplex或gurobi并把输出信息重定向到文件方便回溯调试。4. 系统建模电热综合能源系统的约束与参数怎么定4.1 系统组成与设备运行域算例系统不宜太复杂适合复现和迭代调试。我通常采用改进的IEEE 6节点电力系统耦合一个4节点热力系统。主要设备包括两台常规燃煤机组G1、G2作为第一阶段主力调节电源一台抽凝式CHP机组电热联合运行可行域是典型的梯形区域一台燃气锅炉GB作为热力备用一个储热罐HS承担热负荷转移功能风电场WF输出功率为不确定量电力负荷和热负荷其中热负荷也带有一定随机性。关键设备运行约束要写对尤其是CHP的运行域。抽凝式CHP的可行电出力与热出力关系可以近似表示为[ P_{\text{CHP}} \in [\underline{P}{\text{CHP}}(H{\text{CHP}}),\ \overline{P}{\text{CHP}}(H{\text{CHP}})] ]即电出力上下界都是热出力的线性函数构成一个梯形。这组约束决定了系统的电热耦合强度——CHP在满足热负荷的同时实际上绑架了部分电出力这是电热综合能源系统调度难度的根源。4.2 不确定性数据生成与场景缩减要体现数据驱动首先要有历史数据场景集。实际项目里我们可能只有预测曲线和误差统计信息因此常通过蒙特卡洛或拉丁超立方采样生成N个等概率样本。例如风电出力可表示为预测值加误差[ P_{\text{wind},k} P_{\text{wind,forecast}} \varepsilon_k,\quad \varepsilon_k \sim \mathcal{N}(0, \sigma_w^2) ]热负荷同样处理。要注意样本数量N不能太小——范数约束模糊集的经验分布是以N个场景为支撑的N太小则经验分布可信度低模糊集半径也难以合理化。我一般取N200~500。如果N太大第二阶段LP的求解规模会上升不过现在CPLEX解几百个场景的LP是没问题的影响不大。场景生成好后对每个场景做归一化确保数值稳定性。这是个大坑——不归一化的话热负荷数值可能上千风电出力数值却只有几十两者一起进入模糊集的范数距离计算时风电的偏差会被热负荷的绝对值掩盖导致模糊集形同虚设。4.3 存储与参数标定给出一个实用的参数配置表可以直接用于复现。这个表是我在调参中积累的默认档位适合做一个基准算例后续再根据你的研究需要调整。参数名称数值说明场景数 N200历史数据采样数θ1 (1-范数半径)0.6总概率偏差上限θ∞ (∞-范数半径)0.05单场景概率偏差上限收敛间隙 ε1e-4CCG迭代停止条件弃风惩罚系数50 $/MWh第二阶段目标惩罚切负荷惩罚系数500 $/MWh远高于正常成本防止随意切负荷储热罐容量200 MWh储热罐可存储热量上限CHP最大电出力200 MW—燃气锅炉最大热出力150 MWth—注意θ1和θ∞不是拍脑袋定的需要通过置信水平反推。常用的经验公式是在置信水平α下θ1满足(P{|p-p_0|_1 \le \theta_1} \ge \alpha)。这需要用到概率不等式估算。工程上最简单的方式是换几组θ值跑一遍观察结果成本的灵敏度。θ越大越保守成本越高选在成本急剧增加之前的拐点附近通常是一个合理选择。这一点我会在后面详细展开。5. 算例结果解读怎么证明分布鲁棒比随机规划更好5.1 对比实验设计想让论文或报告有说服力就要做对比。通常设置三组方案随机规划SP直接用经验分布(p_0)作为唯一分布不走模糊集传统鲁棒优化RO不确定量落在区间内不区分概率最坏情况下取极端场景两阶段分布鲁棒DRO本文所用模型模糊集由1-范数和∞-范数约束构造。每组方案跑完后对比三个指标总成本、第一阶段计划、第二阶段最坏期望成本。实际算例跑出来的典型结果是SP的总成本最低但它的风险也最大——一旦真实分布偏离经验分布实际运行成本会急剧上升RO最保守总成本明显偏高第一阶段几乎为所有可能场景都预留了大量备用导致正常日子浪费DRO介于两者之间成本只比SP高几个百分点却能保证在模糊集范围内不出现成本失控。这个结果其实是符合直觉的DRO用数据的概率信息换取了鲁棒性而且通过调节θ1和θ∞可以连续地在SP和RO之间移动——θ越小方案越接近SPθ越大方案越接近RO。这种可调性在工程决策中非常宝贵。5.2 置信度分析为什么不真实分布怎么应对再做一个更有说服力的验证通过蒙特卡洛模拟随机生成大量真实分布不一定等于经验分布然后分别把SP、RO、DRO三套方案放到这些真实分布下计算实际运行成本统计成本分布。结果一般会显示SP方案的成本分布尾部很长——偶尔出现极端高成本RO方案的成本分布集中但均值很高DRO方案的成本分布尾部截断明显优于SP均值又显著低于RO。这种尾部风险控制正是分布鲁棒真正的卖点不是压最低平均成本而是在不严重牺牲经济性的前提下控制最坏情况下的实际损失。如果你的算例结果呈现这个趋势说明你的模型和代码写对了。5.3 参数灵敏度分析θ怎么选才合理θ1和θ∞的选取没有绝对标准但做好灵敏度分析是论文必备内容。我也是在反复实验中发现θ1对结果的影响通常比θ∞更明显因为1-范数约束控制整体概率偏移的总量而∞-范数主要限制单个场景的极端偏移。实际操作时可以固定θ∞不变把θ1从0.2、0.4、0.6、0.8、1.0依次增大观察总成本和迭代次数的变化。典型趋势是θ1增大到某个阈值之前总成本缓慢上升超过某个阈值后成本可能大幅跳增此时模糊集过大、方案过度保守。工程上选择阈值之前的拐点作为推荐值比如0.6。同时θ1增大会让CCG迭代次数增加因为子问题更容易找出不同的极端场景。这个成本-鲁棒性-计算量的综合权衡需要在报告里用图表呈现结论才扎实。6. Matlab代码调试中的常见问题与排查技巧6.1 问题速查表现象可能原因解决办法子问题LP求解失败或不可行模糊集为空θ1太小或θ∞太小导致p无法满足所有约束增大θ1或θ∞检查sum(p)1与范数约束是否冲突CCG迭代次数过多50模糊集半径过大场景数量过多收敛间隙太小适当减小θ1对场景数做缩减放宽ε到1e-3上界持续不下降主问题没有加入正确的割约束第二阶段成本计算与子问题不一致核对主问题新增约束是否使用了子问题返回的p_star检查Q_k表达式是否一致求解结果出现负成本或异常大数目标函数符号写反惩罚系数设置过大导致数值问题逐项检查目标函数调整惩罚系数到合理范围YALMIP报Unable to instantiate solver求解器路径未配置或版本不兼容在yalmip/operators中检查求解器是否支持LP/MILP重新添加CPLEX/Gurobi路径6.2 我调试了几个小时才发现的细节有几个非常不起眼的细节浪费过我大量时间这里专门写出来。第一个是关于YALMIP中repmat的使用。当第二阶段约束里涉及多个场景的表达式时很多人习惯用循环逐个写约束。场景数少还行场景数一多模型构建时间成倍增长。更好的做法是向量化用repmat或直接构造二维变量矩阵一次性写入所有场景的约束。这么做模型构建时间可以从几十秒降到一两秒。第二个是通过value()函数断言检查结果。每次迭代后手动校验一下场景概率p_star是否满足sum(p)1是否都在[0,1]区间内。这个习惯帮我发现过两次模糊集约束写错的问题——一次是因为缺少sum(p)1另一次是因为∞-范数约束误写成了逐元素约束的平方和。第三个是关于求解器日志。CPLEX在求解MILP时日志非常长但关键信息是Best bound和Best integer两个值。建议写一个简单的日志解析脚本自动提取这两个数并绘制收敛曲线能直观判断CCG和MILP的收敛趋势。没有这条曲线你很难判断迭代是在收敛还是卡住了。第四个是一个容易混淆的点1-范数约束中的θ1与∞-范数中的θ∞它们与场景数的关系需要单独校准。当N改变时经验分布p0的各分量p0_k1/N也会变所以同样的θ1和θ∞在不同N下对应不同的模糊程度。做灵敏度分析时务必固定N否则结果不可比。6.3 提升代码运行效率的三个小技巧场景缩减如果初始采样场景数达到1000可以先做k-means聚类或同步回代缩减SBR把场景数降到200~300再进入模糊集模型。这样CCG的迭代次数和子问题规模都会显著降低。热启动Warm StartCCG每次迭代求解主问题时可以以上一轮的解作为初始点YALMIP的assign和optimize的usex0参数可以配合使用。在MILP问题中热启动效果尤其明显。缓存子问题的对偶解在CCG迭代过程中相同场景在相邻迭代中可能被反复求解如果能够缓存每个场景的对偶变量初始值可以加速LP求解。虽然CPLEX自带了一些记忆机制但在自定义实现里显式复用初始点仍然能带来可感知的提升。7. 我个人的几点实操体会这个项目做下来我最大的体会是数据驱动两阶段分布鲁棒优化并不像有些论文写得那么玄乎它本质上就是用数据约束概率分布的不确定性范围然后在这个范围内做最坏情况决策。它的难点反而在细节——模糊集怎么表达、对偶怎么转、CCG的上下界怎么算、YALMIP的向量化怎么写每一步都有不止一个坑。但把这些坑踩平之后你会发现这套框架在电热综合能源系统里特别好用因为CHP的强耦合特性和储热装置的时间转移特性恰恰能让两阶段结构发挥出最大价值。最后再分享一个小技巧当你面对自己的算例结果不知道如何解释时先不要急着调参数而是用蒙特卡洛法把SP、RO、DRO三套方案放到真实分布下做一次批量模拟画一张成本分布箱线图。这张图基本能一眼看出DRO的价值也是论文中最容易打动审稿人的一张图。很多时候方案好坏的差异并不在第一阶段的计划表上而在于它面对真实不确定性时的耐受力——这个道理做调度的人应该都懂。
