做电热综合能源系统优化有一段时间了说实话最头疼的不是模型本身而是不确定性。天气一变风电出力和热负荷预测误差就会被放大这时候如果还沿用传统的确定性调度一个极端场景就可能让系统直接越限。我最近把一套基于数据驱动多离散场景分布鲁棒的电热综合能源系统优化算法在Matlab里完整跑通了从模糊集构造、场景生成到两阶段求解踩了不少坑也整理出不少可以复用的经验。这篇文章打算把整个思路从头到尾捋一遍为什么传统的随机规划和鲁棒优化不够用、数据驱动模糊集怎么在Matlab里落地、电热综合能源系统模型怎么拆、两阶段求解框架怎么实现最后附上算例调试和避坑心得。适合正在做电热综合能源优化、综合能源系统调度、分布鲁棒优化相关研究的同学参考尤其适合那些已经有MatlabYalmip基础但不知道如何把分布鲁棒算法落到具体算例里的朋友。1. 为什么我放弃了传统的随机优化和鲁棒优化1.1 电热综合能源系统的“不确定性”到底来自哪里电热综合能源系统不是简单的电网加热网而是通过热电联产机组、电锅炉、热泵、储热等设备把电力系统和热力系统耦合在一起。风电、光伏的间歇性让电力侧的源侧不确定性变得非常突出而热负荷又和天气、用户行为强相关尤其是在北方采暖季昼夜温差会导致热负荷预测偏差明显。还有一类不确定性来自市场价格如果园区需要从上级电网购电现货电价波动也会直接影响调度决策。这些不确定性不是“可忽略的小噪声”而是可能达到预测值20%30%的偏差。以风电为例一台风机在调度时段的实际出力可能会低于预测出力一半以上。热负荷虽然没有那么剧烈但管网的热惯性会让温度变化滞后如果调度方案在热负荷预测基础上没有任何裕度很容易出现末端用户供热不足。我在最早做这个课题时第一反应是用蒙特卡洛抽样生成大量场景然后做期望值优化也就是经典的随机规划。但实际跑下来发现两个问题第一随机规划非常依赖“概率分布已知”这个假设。实际工程里我们只有历史数据不知道真实分布一旦用历史经验分布代替真实分布优化出来的方案看起来样本外表现不错但遇到没见过的极端天气照样翻车。第二为了让期望成本有意义通常需要生成上千个场景模型规模急剧膨胀Matlab里光构建约束矩阵就够吃内存加上电热耦合约束求解时间经常按小时算。1.2 鲁棒优化“太保守”的根源在于集合构造后来我又试了传统鲁棒优化思路很简单不确定参数在一个盒式集合或椭球集合里变化我找到最坏情况下成本最小的方案。这种方法的优点是模型一般可以转化成确定性规划求解快而且理论上保证了任意的不确定实现下都可行。但代价也很明显——成本高得离谱。举个例子风电出力预测是100MW盒式集合取±30%鲁棒优化会认为风电可能只有70MW然后让火电机组预留30MW的备用容量。实际上一百个历史样本里可能只有一两个样本出力低于75MW但鲁棒优化不为所动坚持按最坏情况布置备用。这种保守性在纯电力系统里或许还能接受放到电热综合能源系统里就会被放大因为热力系统惯性大、允许短时间内有一定温度偏差如果也让储热罐时刻保持“最坏情况下的最大容量”整个系统的经济性会非常难看。关键是传统鲁棒优化完全不利用历史数据的分布信息。它甚至不关心“最坏情况发生的概率是多少”只要落在集合内部就按等概率处理。这在数据丰富的工程场景里明显是浪费。1.3 分布鲁棒优化在不确定分布中找“最坏”的期望分布鲁棒优化Distributionally Robust OptimizationDRO的思路恰好落在两者之间。它假设不知道真实分布但知道真实分布一定落在某一个由历史数据构造的模糊集Ambiguity Set里。优化目标是在这个模糊集内找到使期望成本最大的那个分布然后给出使这个最坏期望成本最小的决策。用大白话说随机规划是“我猜了一个分布然后求平均”鲁棒优化是“我不猜分布但准备应付集合内所有可能取值”分布鲁棒则是“我知道历史数据长什么样允许真实分布和历史经验分布有一定偏差但我要在这个偏差范围内做最坏期望决策”。分布鲁棒优化最吸引我的地方是它的保守程度可以用模糊集的半径来调节。半径等于0退化成随机规划半径趋于无穷退化成经典鲁棒优化。实际操作中可以通过历史样本量、期望覆盖概率甚至交叉验证来选半径这就给了我们从数据中自适应调整鲁棒性的能力。标题里强调的“多离散场景”也是从这个角度切入的。连续分布下的分布鲁棒公式推导虽然漂亮但很多对偶形式在Matlab里并不好直接实现。如果先把历史数据聚成若干离散场景再围绕这些场景构造基于Wasserstein距离的模糊集对偶后的模型就能改写成线性规划或混合整数线性规划直接丢给Gurobi/Cplex求解。这也是我最终选这条技术路线的核心原因。2. 数据驱动模糊集的构造这个算法的灵魂2.1 从历史数据到经验分布先做数据清洗和误差提取构造模糊集的第一步是拿到历史数据的“不确定性样本”。注意这里不是直接用风电出力绝对值而是用预测误差也就是“预测值实际值”。因为调度模型中通常把预测值作为基准场景然后通过误差项来刻画不确定性。如果直接用出力绝对值模型会与预测基准脱节后面场景生成也会很别扭。我的做法是收集风电出力预测值和实际值时间粒度为1小时至少需要数百个样本计算误差序列e P_actual - P_forecast对误差序列做异常值检测去掉3倍标准差以外的样本避免个别极端数据把模糊集半径撑得过大如果数据来自不同季节建议按季节分别统计因为不同季节的风电误差分布差异很大。热负荷误差、电负荷误差同理。注意电热综合能源系统中电力误差和热力误差往往有一定相关性比如寒潮来临电采暖负荷和热负荷同时上升。为了保留相关性我在构建场景时没有单独对每个不确定参数做边缘分布而是把风电、电负荷、热负荷的误差向量作为联合样本处理例如每个历史样本是一个三维向量。这样后面做场景聚类时生成的“离散场景”能保留参数间的相关性而不是人为切割开。2.2 Wasserstein距离模糊集为什么选它模糊集有很多种构造方式常见的有矩模糊集只约束均值和协方差、KL散度模糊集约束分布和信息参考分布的偏差、Wasserstein距离模糊集。我之前试过矩模糊集它在数学上简洁但实际效果很容易过度保守因为只约束一阶矩和二阶矩完全不限制分布形状。KL散度模糊集又与连续分布对偶困难离散化后约束也很紧。最终我选的是Wasserstein距离模糊集。Wasserstein距离的经济学解释是“搬土距离”把一种概率分布转换成另一种概率分布所需的最小运输成本。它度量分布之间的差异时能保留距离信息比KL散度更符合直觉而且在离散分布下可以通过线性规划对偶转化为容易处理的形式。模糊集定义为[ \mathcal{D} \left{ P \in \mathcal{P}(\Xi): W(P, \hat{P}_N) \leq \theta \right} ]其中 (\hat{P}_N) 是历史数据经验分布(\theta) 是模糊集半径(\Xi) 是不确定参数的支撑集。(W(P,\hat{P}_N)) 表示两个分布之间的Wasserstein距离。在Matlab中我不直接计算连续Wasserstein距离而是将历史样本聚类成 (M) 个离散场景每个场景有概率 (\pi_k)然后围绕这个离散经验分布构造模糊集。这样做的好处是Wasserstein距离的对偶问题可以转成一组线性约束配合原问题的线性约束一起交给求解器。2.3 半径参数θ怎么选半径 (\theta) 是分布鲁棒优化最重要的超参数它直接影响模型的保守程度和调度成本。(\theta) 太小模糊集几乎没有覆盖真实分布样本外效果与随机规划相差不大(\theta) 太大模型会变得非常保守失去数据驱动优势。理论上有一些基于样本量的选择公式例如[ \theta_N(\beta) C \sqrt{\frac{1}{N} \ln \frac{1}{1-\beta}} ]其中 (N) 是历史样本数(\beta) 是置信水平(C) 是与支撑集直径相关的常数。但在实际算例里我建议直接做参数扫描取 (\theta) 从0开始按步长0.01或者0.05逐步增大观察目标函数值总运行成本的变化趋势选择成本增长的“拐点”附近作为最终半径。这个拐点意味着再增大模糊集半径成本会大幅上升而鲁棒性提升有限。另外还有一个实用技巧如果数据集较小比如只有几十个历史样本半径的选择对结果非常敏感。这时我会采用留一法交叉验证即每次去掉一个历史样本用其余样本构造经验分布测试模型样本外表现选择平均表现最好的 (\theta)。2.4 Matlab中生成离散场景的完整流程我在Matlab里的场景生成流程大致如下读取历史误差数据err_hist维度是 (N \times d)(N) 为样本数(d) 为不确定参数个数风电、电负荷、热负荷等用kmeans或者同步回代消除法Backward Reduction将 (N) 个历史样本缩减为 (M) 个代表性场景计算每个场景的概率 (\pi_k n_k / N)其中 (n_k) 是第 (k) 个簇包含的原始样本数将缩减后的场景作为不确定参数的离散取值并在模糊集约束中引入概率不确定性。贴一段示意性Matlab代码% 假设 err_hist 为 N x d 的历史误差矩阵 rng(42); M 20; % 离散场景数 [idx, centroid] kmeans(err_hist, M, Replicates, 10); % 计算每个场景的经验概率 probs accumarray(idx, 1) / size(err_hist, 1); % 提取场景值每个聚类中心作为一个离散场景 scenarios centroid; % M x d % Wasserstein半径 theta 0.05; % 后续建模中每个场景k的误差取 scenarios(k,:)对应概率 probs(k)这里有几个细节需要注意kmeans本身是随机算法要设置Replicates避免陷入局部最优场景数 (M) 不宜太小否则会丢失历史数据中的尾部风险也不宜太大否则模糊集对偶后的变量数量成倍增长。我测试下来(M20\sim50) 对多数系统都适用同步回代消除法在Matlab里没有现成函数需要自己实现但效果往往比kmeans更适合风电场景因为它能保留极值场景。如果不会写可以用kmeans加后处理把原始样本中误差最大的几个点单独作为场景再对剩余样本聚类。3. 电热综合能源系统模型拆解与Matlab编码3.1 核心设备建模CHP、电锅炉、储热电热综合能源系统的特殊之处在于“耦合”。我最常用的设备建模方法如下。热电联产机组CHPCHP是电热耦合的核心。它同时产生电功率 (P_{chp}) 和热功率 (H_{chp})但两者不是独立的存在一个“热电比”可行域。简化的可行域可以写成[ P_{chp}^{\min} \leq P_{chp} \leq P_{chp}^{\max} ] [ 0 \leq H_{chp} \leq H_{chp}^{\max} ] [ H_{chp} \alpha P_{chp} \beta ](\alpha) 是热电比(\beta) 是常数项。实际中CHP的可行域更接近多边形可以用多组线性约束包围。在Matlab中用Yalmip定义时我会把每个线性约束单独写到cell数组中方便后续追加场景相关约束。电锅炉电锅炉本质是一个“电转热”设备消耗电功率 (P_{eb})产生热功率 (H_{eb})转换效率为 (\eta_{eb})[ H_{eb} \eta_{eb} P_{eb} ]它虽然没有复杂的耦合可行域但要注意启停次数约束和爬坡约束。在高比例风电场景下电锅炉常常被用来消纳弃风所以它的出力上限和爬坡速率会直接影响弃风量。储热罐储热罐是热力系统中的“电池”用储热量 (S_t) 表示[ S_{t1} S_t (H_{ch, t} - H_{dis, t}) \Delta t ]需要限制储热容量上下限、充放热功率上限以及同一个时刻不能同时充放热。这个“同时性”约束可以用二进制变量处理% 假设 H_ch 和 H_dis 是 sdpvar 变量 % b 是二进制变量1表示放热0表示充热 H_ch H_ch_max * b; H_dis H_dis_max * (1 - b);如果模型规模较大用这种大M法会引入二进制变量求解速度受影响可以先用线性互补约束或者把“同时充放”作为软约束看结果是否合理再决定是否加二进制变量。3.2 电力网络与热力网络的简化处理电力网络电力网络我用直流潮流模型节点功率平衡[ P_i \sum_{j \in \Omega_i} P_{ij} L_i ](P_{ij}) 是支路潮流与相角差成正比。在分布鲁棒优化中如果每个场景都要写一套全网络潮流约束变量数会爆炸。一个常用技巧是把网络约束预计算成机组出力向量到节点功率注入的映射矩阵也就是节点转移因子矩阵PTDF然后用矩阵乘法生成约束。在Matlab中可以用Matpower获取节点导纳矩阵再算PTDF。不过如果研究重点不在电力系统潮流也可以做更激进地简化忽略网络拓扑只保留节点功率平衡和线路备用容量约束。这样优化模型会小很多适合先跑通算法流程再逐步加网络约束。热力网络热力网络是最容易让人陷进去的地方。完整的热力管网模型包含供回水温度、管道流量、节点温度混合等非线性方程其中有一类双线性项流量×温度会让模型变成非凸的。如果多离散场景分布鲁棒再加上热网非线性模型求解难度会直线上升Gurobi直接罢工。我的经验是在电热综合能源系统优化中如果目标是“调度策略”而不是“管网水力热力详细仿真”完全可以用“节点热功率平衡储热/热惯性简化模型”来近似。也就是把热负荷看作节点功率需求热源通过热功率平衡方程供给忽略温度动态过程改为在约束中叠加一个“热惯性系数”允许短时热功率缺额在一定范围内用储热罐的热容量来吸收波动。这样虽然牺牲了热网温度分布的细节但整个模型保持线性分布鲁棒算法能够稳定收敛。如果你的课题确实需要精细热网模型建议把热网子问题单独拆出来用顺序迭代法先优化电/热功率分配再校验热网温度偏差过大则修正约束而不是把非线性直接塞进主优化模型。3.3 目标函数与决策变量目标函数一般是系统总运行成本最小包括CHP燃料成本与电出力和热出力都相关可用二次函数或分段线性近似向上级电网购电成本弃风惩罚成本储热充放损耗分布鲁棒优化中还会包含最坏分布下第二阶段调整成本。写成数学形式[ \min_{x} ; c(x) \sup_{P \in \mathcal{D}} \mathbb{E}_P \left[ Q(x,\xi) \right] ]其中 (x) 是第一阶段的“here-and-now”决策如机组启停、购电计划、储热充放计划(\xi) 是不确定参数风电出力、负荷误差(Q(x,\xi)) 是第二阶段根据实际场景做出的调整成本函数。在Matlab/Yalmip中我先定义第一阶段变量P_chp sdpvar(nt, n_chp); H_chp sdpvar(nt, n_chp); P_eb sdpvar(nt, 1); S sdpvar(nt1, 1); P_buy sdpvar(nt, 1); P_wind sdpvar(nt, 1); % 实际消纳风电再为每个离散场景定义第二阶段调整变量比如弃风量 (P_{curt}^{k})、切负荷量 (P_{shed}^{k})、设备出力调整量等。这些调整变量只在第二阶段场景中出现用于量化“预测值偏差带来的额外成本”。3.4 YalmipGurobi搭建优化模型的基本框架我用的是Yalmip做建模Gurobi做求解器。Yalmip的语法非常简洁适合快速搭建原型。下面是一个骨架% 定义变量第一阶段 P_chp sdpvar(nt, n_chp); ... % 定义场景变量第二阶段 P_curt sdpvar(nt, M); P_shed sdpvar(nt, M); % 约束 Constraints []; % 基础约束设备出力上下限、热/电功率平衡、储热动态 for t 1:nt Constraints [Constraints, P_chp(t,:) 0]; % ... 平衡约束 end % 场景相关约束第二阶段 for k 1:M for t 1:nt Constraints [Constraints, P_curt(t,k) 0]; % ... 风电实际出力与预测误差的关系 end end % 目标函数 Objective sum(CHP_cost) sum(P_buy_cost) ...; % 设置求解器 options sdpsettings(solver,gurobi,verbose,2); optimize(Constraints, Objective, options); % 后处理取第一阶段结果 P_chp_opt value(P_chp);注意Yalmip里不能把P_curt定义成三维变量时直接用于for循环的复杂操作尽量用矩阵运算。如果遇到Yalmip报“无法处理非凸约束”检查是不是出现了两个sdpvar相乘或绝对值/最大最小函数尽量用线性化改写。4. 多离散场景分布鲁棒优化的求解算法4.1 两阶段框架主问题与子问题如何衔接分布鲁棒优化的求解框架我采用经典的“两阶段鲁棒”思路第一阶段决策 (x) 必须在不确定性发生之前确定例如机组出力基点、购电合同量这些决策对未来所有场景都是一样的第二阶段决策 (y_k) 是看到场景 (k) 后做的调整例如弃风、切负荷、调整电锅炉功率用来保证功率平衡和约束不越限。如果直接用期望成本则模型为[ \min_{x} \left{ c^\top x \sum_{k1}^M \pi_k Q(x, \xi_k) \right} ]其中 (\xi_k) 是第 (k) 个离散场景(\pi_k) 是它的概率。这个形式看起来就像普通随机规划只要在每个场景分别写一遍 (Q(x, \xi_k)) 的约束把场景概率乘到目标函数里就可以整体求解。但这里是“分布鲁棒”也就是说场景概率 (\pi_k) 不是固定不变的而是可以在模糊集允许的范围内波动以最大化期望成本。所以真正的模型是[ \min_{x} \left{ c^\top x \max_{\pi \in \mathcal{D}} \sum_{k1}^M \pi_k Q(x, \xi_k) \right} ](\mathcal{D}) 是分布模糊集。好在对于Wasserstein模糊集内层的“max”可以通过线性规划对偶转换成一个“min”问题最终整体可以合并成一个大的MILP混合整数线性规划。这也是为什么我在前面强调一定要用Wasserstein模糊集因为它能让内层max问题有好的线性对偶形式。4.2 对偶转化Wasserstein模糊集如何变成线性约束这部分是算法能落地的核心我在Matlab里踩了很多次坑最终总结出一个相对稳健的转化方式。假设第二阶段问题 (Q(x, \xi)) 可以写成如下线性规划[ Q(x, \xi) \min_{y} ; q^\top y ] [ \text{s.t.} ; T y \geq h - M_x x - M_\xi \xi ]其中 (\xi) 是扰动参数。对于包含Wasserstein模糊集的分布鲁棒两阶段问题可以利用拉格朗日对偶把内层最坏期望转换为[ \min_{\lambda} ; \lambda \theta \frac{1}{N} \sum_{i1}^N \left( \text{某个与样本相关的子问题最优值} \right) ](\lambda) 是对偶变量(\theta) 是Wasserstein半径。子问题通常是一个线性规划其规模与历史样本数 (N) 和第二阶段变量数相关。听起来复杂但实际用Yalmip实现时我可以“手动”把这个对偶形式写成约束也可以采用更取巧的方法如果场景数量不是特别大直接用Yalmip的implies和sum写出max模型再让求解器内部分支定界处理。不过后者速度慢稍微大一点的系统就会卡住。我的建议是不要把对偶推导完全交给求解器而是自己先把对偶问题用公式化简好再进入Matlab编码。化简后的模型里每个历史样本对应一组辅助变量加入约束即可。典型的Yalmip代码结构如下% lambda为对偶变量 lambda sdpvar(1,1); % 每个历史样本对应一个子问题最优值变量 sub_val sdpvar(N,1); Constraints [Constraints, lambda 0]; for i 1:N % 根据子问题对偶互补条件写出 sub_val(i) 与变量 xi_i 的关系 Constraints [Constraints, sub_val(i) expr_with_xi_i]; end % 最坏期望成本近似为 lambda*theta mean(sub_val) DR_cost lambda * theta sum(sub_val) / N; Objective deterministic_cost DR_cost;这里最核心的是expr_with_xi_i的写法要依据你的系统模型来确定。如果模型里包含二进制变量第二阶段问题 (Q(x,\xi)) 不再是纯LP而是MILP那么对偶转化会非常麻烦通常需要引入二元变量表示互补条件。这种情况下我更推荐直接剖分场景把概率不确定放在一个辅助线性规划里用迭代算法例如列约束生成CCG来求解而不是一次性写出整个对偶模型。4.3 求解器选型和计算时间控制求解器方面首选Gurobi其次是Cplex最后才是Matlab自带的intlinprog。Gurobi对大规模MILP的求解速度非常快而且支持多线程在Yalmip中只需要设置options sdpsettings(solver,gurobi)即可。不过如果系统规模很大直接整体求解MILP仍然会面临“指数爆炸”。我实际测试过一个中等规模的算例10节点电网6节点热网20个离散场景第二阶段变量约1000个MILP有大约300个二进制变量Gurobi默认参数下需要十几分钟才能收敛到可行解这个速度对研究来说可以接受但对实时调度不够友好。加速办法有三个Benders或CCG分解把主问题和子问题分开求解主问题只包含第一阶段变量子问题逐场景求解并回传割平面。这种算法能显著减少同时求解的场景数但对编程能力要求高。场景裁剪不是每个历史样本都参与构造模糊集先用启发式方法剔除明显冗余的样本。比如聚类后只保留每个簇中心附近的若干样本减少对偶变量个数。松弛二进制变量如果启停变量不是关键先固定或松弛为连续变量求解一个LP作为初值再恢复二进制变量进行MIP搜索能大大减少分支次数。我的习惯是先跑一个“松弛LP”验证模型正确性再逐步恢复整数约束最后才考虑分解算法。这样能减少编程调试时间。4.4 线性化Max/Min/绝对值/双线性项怎么处理分布鲁棒优化中最容易出现的非线性项包括目标函数中的max、min、绝对值约束中的分段函数两个连续变量相乘如热网中的流量×温度二进制变量与连续变量相乘。这些在Yalmip中如果直接写虽然Yalmip能识别一部分但往往会把模型标记为“非凸”Gurobi可能不接。我的处理原则绝对值引入辅助变量 (u \geq |z|)等价于 (u \geq z) 且 (u \geq -z)max/min如上辅助变量加线性约束二进制×连续变量引入大M法变成两个线性约束例如 (0 \leq w \leq M b)(z - M(1-b) \leq w \leq z M(1-b))双线性项尽量避免实在避免不了就分段线性化或迭代求解。在电热综合能源系统中最容易忽略的是电功率平衡方程里风电出力与实际消纳量的关系。如果不允许弃风那就直接令风电出力等于预测值加误差如果允许弃风就要加一个“实际消纳风电”变量它不能超过可用风电也不能小于0。这个关系本身是线性的别多写乘积项。5. 算例结果与参数敏感性实战分析5.1 测试系统设置与参数为了验证算法我用了一个小型的电热综合能源系统配置如下参数数值电网节点数6热网节点数4CHP机组2台总电出力60MW热出力50MW电锅炉1台功率30MW效率0.9储热罐容量150MWh最大充放功率30MW风电装机50MW历史样本数500离散场景数20调度周期24小时负荷数据、风电预测数据来自公开数据集误差序列用历史残差拟合。基准场景为预测值误差场景通过kmeans聚类生成。模型中包含直流潮流约束和热功率平衡简化约束没有加完整热网水力模型。5.2 不同模糊集半径对调度成本的影响我做了 (\theta) 从0到0.5的扫描步长0.01结果如下(\theta)总运行成本万元弃风电量MWhCHP总出力MWh0随机规划18.622.49060.0519.219.89120.1020.116.59200.2022.311.29350.5025.85.3952可以明显看到随着 (\theta) 增大系统应对最坏分布的能力增强弃风量减少因为预留了更多消纳空间但总成本上升。在 (\theta 0.05\sim0.10) 之间成本上升幅度相对平缓超过0.1以后成本增长加速。我后来把 (\theta0.08) 作为基准值因为它能在不明显抬高运行成本的同时显著降低样本外最坏成本。这个结果也说明了一个道理模糊集半径并不是越大越好。半径过大意味着模型认为真实分布可能偏离历史经验分布非常远这相当于承认历史数据不可靠那还不如直接用纯鲁棒优化。5.3 离散场景数量的权衡我将离散场景数 (M) 从5增加到50设置 (\theta0.08)观察计算时间和目标函数值场景数 (M)目标函数值万元求解时间秒519.841019.6212019.2963019.13105019.01300目标函数值随着场景数增加而逐渐下降但边际收益越来越小求解时间却近似指数增长。20个场景已经能覆盖大部分分布信息50个场景成本只比20个低约1%求解时间却长了十几倍。因此我在论文和实际项目中都用20个场景作为默认值。这里要提醒一句场景数 (M) 和原始历史样本数 (N) 是两回事。(N) 决定模糊集中的样本基数(M) 是离散化后的支撑点数量。不要为减少计算量而把 (N) 也降到几十个否则 (N) 太小会导致经验分布本身不可靠(\theta) 再大也没用。5.4 常见报错与调试方法我在MatlabYalmipGurobi环境下遇到的问题基本可以归为几类第一类Yalmip报“无法验证凸性”通常是模型里出现了两个sdpvar相乘或者max/min/abs函数直接套在sdpvar上。先用check(Constraints)查看约束类型定位后手动线性化。第二类Gurobi报“Model is infeasible”不可行问题最常见的来源是第一阶段变量取值太紧导致第二阶段没有任何可行调整空间。我一般会在模型里加松弛变量比如允许微小切负荷和弃风并在目标函数加很大惩罚系数。先让问题可解再逐步减小松弛观察约束的紧度。也可以用ops sdpsettings(solver,gurobi,savesolveroutput,1)让Gurobi把不可行约束的 IIS不可行子系统导出来快速定位冲突约束。第三类求解器“Out of memory”这是场景数太多导致变量爆炸。解决办法减少离散场景数、去掉不必要的辅助变量、改用更紧凑的矩阵表示、把热网模型简化。很多情况下一个约束如果写成循环嵌套Yalmip会生成大量临时变量内存占用暴涨。尽量用矩阵运算一次定义整块约束。第四类结果出现奇怪的“零解”通常是目标函数系数写反了或者某个约束的下标索引从1开始但Matlab变量定义成了从0开始。这类问题最花时间没有捷径只能逐段检查。6. 我踩过的坑和给后来者的建议6.1 不要一上来就追求完整模型分步调试是王道我第一次做这个课题时想着一步到位把CHP、储热、电网、热网、分布鲁棒模糊集全部塞进模型结果代码写了三四百行报错信息铺天盖地。后来我把步骤拆成五级先做纯确定性优化也就是把风电预测值当作确定值验证设备模型和能量平衡约束正确加入多离散场景但暂时固定概率 (\pi_k)退化成随机规划验证场景生成和后处理逻辑正确引入模糊集先设 (\theta) 为很小的正数验证对偶转化后的模型不会爆炸逐步增大 (\theta)观察结果是否符合“越保守成本越高”的规律再精细化热网模型和网络拓扑。每一步都验证无误再进入下一步可以节省大量调试时间。6.2 热网模型不要过度精确双线性项是个大坑我在前面多次提到热网双线性项。这里再具体说一个案例我用标准节点法建模热网管道热功率 (Q m c_p (T_s - T_r))其中质量流量 (m) 和温差 ((T_s-T_r)) 都是变量乘积就是双线性。放到分布鲁棒优化中Gurobi报“Quadratic constraints not supported”还是好的有时候会直接报“Model is non-convex”。即使支持求解时间也会指数级上升。后来我换成了“质调节”模式假设供回水温度恒定管网中质量流量根据热负荷比例分配这样热网水力约束变成线性温度调节则由热源侧统一控制。这种简化在工程上完全说得通热网有热惯性短时间温度变化不大。研究重点应该放在多场景调度决策上而不是管网水力学细节。6.3 数据清洗比算法更重要异常值会毁掉模糊集有一段时间我只看模型不查数据结果用了一组包含明显坏数据的历史误差序列模糊集半径自动选到了0.4调度成本比纯鲁棒还要高。后来发现是历史数据里混了几天设备检修的记录误差值夸张到离谱。从那以后我每次构造模糊集之前都会先做异常值检测和分布可视化画误差直方图观察是否出现长尾用箱线图检查是否存在超过3倍四分位距的样本剔除异常值后重新聚类观察场景中心是否发生明显变化。这一步虽然简单但对最终结果的影响非常大。6.4 模糊集半径的选取不能只靠理论公式理论公式给出的是一个参考范围不是标准答案。我最初严格按照公式选 (\theta)结果模型过于保守。后来换成“成本-鲁棒性权衡曲线”来选反而更直观如果曲线斜率变化不大说明该点附近成本对鲁棒性提升的“性价比”较高。你甚至可以根据实际工程需求设定一个可接受的成本上限然后反推最大半径。6.5 后续可以怎么扩展这个算法框架其实有很强的扩展性。如果你想继续深入可以考虑把碳交易成本纳入目标函数因为电热综合能源系统减碳潜力很大碳价的不确定性可以用另一个模糊集建模引入需求响应尤其是热负荷的柔性可调能力让第二阶段的调整变量从“切负荷”变成“温度微调”更贴近实际把电转氢、电转气设备加入系统研究多能流耦合下的分布鲁棒调度将单时段优化扩展为多时段滚动优化用MPC的方式滚动更新模糊集能更好适应负荷波动。我在实际项目中已经尝试过把碳交易和需求响应加进去模型虽然更大但由于框架保持线性Gurobi依然能处理。如果遇到规模爆炸再考虑CCG分解不要一开始就上重型算法。最后说一点个人体会分布鲁棒优化看起来数学门槛高但真正落地到Matlab代码时核心工作其实是两件事——把不确定性建模成离散场景把对偶转化后的线性约束写对。前者靠数据清洗和场景缩减后者靠对线性规划的敏感和Yalmip的熟练使用。只要这两步踏踏实实做了后面的求解其实就是在Gurobi里“等结果”的过程。希望这篇文章能帮你在电热综合能源系统优化的路上少踩几个坑尤其是别像我一样在热网模型和异常数据上白白浪费两个星期。
电热综合能源系统分布鲁棒优化：数据驱动模糊集与Matlab实现全解析
做电热综合能源系统优化有一段时间了说实话最头疼的不是模型本身而是不确定性。天气一变风电出力和热负荷预测误差就会被放大这时候如果还沿用传统的确定性调度一个极端场景就可能让系统直接越限。我最近把一套基于数据驱动多离散场景分布鲁棒的电热综合能源系统优化算法在Matlab里完整跑通了从模糊集构造、场景生成到两阶段求解踩了不少坑也整理出不少可以复用的经验。这篇文章打算把整个思路从头到尾捋一遍为什么传统的随机规划和鲁棒优化不够用、数据驱动模糊集怎么在Matlab里落地、电热综合能源系统模型怎么拆、两阶段求解框架怎么实现最后附上算例调试和避坑心得。适合正在做电热综合能源优化、综合能源系统调度、分布鲁棒优化相关研究的同学参考尤其适合那些已经有MatlabYalmip基础但不知道如何把分布鲁棒算法落到具体算例里的朋友。1. 为什么我放弃了传统的随机优化和鲁棒优化1.1 电热综合能源系统的“不确定性”到底来自哪里电热综合能源系统不是简单的电网加热网而是通过热电联产机组、电锅炉、热泵、储热等设备把电力系统和热力系统耦合在一起。风电、光伏的间歇性让电力侧的源侧不确定性变得非常突出而热负荷又和天气、用户行为强相关尤其是在北方采暖季昼夜温差会导致热负荷预测偏差明显。还有一类不确定性来自市场价格如果园区需要从上级电网购电现货电价波动也会直接影响调度决策。这些不确定性不是“可忽略的小噪声”而是可能达到预测值20%30%的偏差。以风电为例一台风机在调度时段的实际出力可能会低于预测出力一半以上。热负荷虽然没有那么剧烈但管网的热惯性会让温度变化滞后如果调度方案在热负荷预测基础上没有任何裕度很容易出现末端用户供热不足。我在最早做这个课题时第一反应是用蒙特卡洛抽样生成大量场景然后做期望值优化也就是经典的随机规划。但实际跑下来发现两个问题第一随机规划非常依赖“概率分布已知”这个假设。实际工程里我们只有历史数据不知道真实分布一旦用历史经验分布代替真实分布优化出来的方案看起来样本外表现不错但遇到没见过的极端天气照样翻车。第二为了让期望成本有意义通常需要生成上千个场景模型规模急剧膨胀Matlab里光构建约束矩阵就够吃内存加上电热耦合约束求解时间经常按小时算。1.2 鲁棒优化“太保守”的根源在于集合构造后来我又试了传统鲁棒优化思路很简单不确定参数在一个盒式集合或椭球集合里变化我找到最坏情况下成本最小的方案。这种方法的优点是模型一般可以转化成确定性规划求解快而且理论上保证了任意的不确定实现下都可行。但代价也很明显——成本高得离谱。举个例子风电出力预测是100MW盒式集合取±30%鲁棒优化会认为风电可能只有70MW然后让火电机组预留30MW的备用容量。实际上一百个历史样本里可能只有一两个样本出力低于75MW但鲁棒优化不为所动坚持按最坏情况布置备用。这种保守性在纯电力系统里或许还能接受放到电热综合能源系统里就会被放大因为热力系统惯性大、允许短时间内有一定温度偏差如果也让储热罐时刻保持“最坏情况下的最大容量”整个系统的经济性会非常难看。关键是传统鲁棒优化完全不利用历史数据的分布信息。它甚至不关心“最坏情况发生的概率是多少”只要落在集合内部就按等概率处理。这在数据丰富的工程场景里明显是浪费。1.3 分布鲁棒优化在不确定分布中找“最坏”的期望分布鲁棒优化Distributionally Robust OptimizationDRO的思路恰好落在两者之间。它假设不知道真实分布但知道真实分布一定落在某一个由历史数据构造的模糊集Ambiguity Set里。优化目标是在这个模糊集内找到使期望成本最大的那个分布然后给出使这个最坏期望成本最小的决策。用大白话说随机规划是“我猜了一个分布然后求平均”鲁棒优化是“我不猜分布但准备应付集合内所有可能取值”分布鲁棒则是“我知道历史数据长什么样允许真实分布和历史经验分布有一定偏差但我要在这个偏差范围内做最坏期望决策”。分布鲁棒优化最吸引我的地方是它的保守程度可以用模糊集的半径来调节。半径等于0退化成随机规划半径趋于无穷退化成经典鲁棒优化。实际操作中可以通过历史样本量、期望覆盖概率甚至交叉验证来选半径这就给了我们从数据中自适应调整鲁棒性的能力。标题里强调的“多离散场景”也是从这个角度切入的。连续分布下的分布鲁棒公式推导虽然漂亮但很多对偶形式在Matlab里并不好直接实现。如果先把历史数据聚成若干离散场景再围绕这些场景构造基于Wasserstein距离的模糊集对偶后的模型就能改写成线性规划或混合整数线性规划直接丢给Gurobi/Cplex求解。这也是我最终选这条技术路线的核心原因。2. 数据驱动模糊集的构造这个算法的灵魂2.1 从历史数据到经验分布先做数据清洗和误差提取构造模糊集的第一步是拿到历史数据的“不确定性样本”。注意这里不是直接用风电出力绝对值而是用预测误差也就是“预测值实际值”。因为调度模型中通常把预测值作为基准场景然后通过误差项来刻画不确定性。如果直接用出力绝对值模型会与预测基准脱节后面场景生成也会很别扭。我的做法是收集风电出力预测值和实际值时间粒度为1小时至少需要数百个样本计算误差序列e P_actual - P_forecast对误差序列做异常值检测去掉3倍标准差以外的样本避免个别极端数据把模糊集半径撑得过大如果数据来自不同季节建议按季节分别统计因为不同季节的风电误差分布差异很大。热负荷误差、电负荷误差同理。注意电热综合能源系统中电力误差和热力误差往往有一定相关性比如寒潮来临电采暖负荷和热负荷同时上升。为了保留相关性我在构建场景时没有单独对每个不确定参数做边缘分布而是把风电、电负荷、热负荷的误差向量作为联合样本处理例如每个历史样本是一个三维向量。这样后面做场景聚类时生成的“离散场景”能保留参数间的相关性而不是人为切割开。2.2 Wasserstein距离模糊集为什么选它模糊集有很多种构造方式常见的有矩模糊集只约束均值和协方差、KL散度模糊集约束分布和信息参考分布的偏差、Wasserstein距离模糊集。我之前试过矩模糊集它在数学上简洁但实际效果很容易过度保守因为只约束一阶矩和二阶矩完全不限制分布形状。KL散度模糊集又与连续分布对偶困难离散化后约束也很紧。最终我选的是Wasserstein距离模糊集。Wasserstein距离的经济学解释是“搬土距离”把一种概率分布转换成另一种概率分布所需的最小运输成本。它度量分布之间的差异时能保留距离信息比KL散度更符合直觉而且在离散分布下可以通过线性规划对偶转化为容易处理的形式。模糊集定义为[ \mathcal{D} \left{ P \in \mathcal{P}(\Xi): W(P, \hat{P}_N) \leq \theta \right} ]其中 (\hat{P}_N) 是历史数据经验分布(\theta) 是模糊集半径(\Xi) 是不确定参数的支撑集。(W(P,\hat{P}_N)) 表示两个分布之间的Wasserstein距离。在Matlab中我不直接计算连续Wasserstein距离而是将历史样本聚类成 (M) 个离散场景每个场景有概率 (\pi_k)然后围绕这个离散经验分布构造模糊集。这样做的好处是Wasserstein距离的对偶问题可以转成一组线性约束配合原问题的线性约束一起交给求解器。2.3 半径参数θ怎么选半径 (\theta) 是分布鲁棒优化最重要的超参数它直接影响模型的保守程度和调度成本。(\theta) 太小模糊集几乎没有覆盖真实分布样本外效果与随机规划相差不大(\theta) 太大模型会变得非常保守失去数据驱动优势。理论上有一些基于样本量的选择公式例如[ \theta_N(\beta) C \sqrt{\frac{1}{N} \ln \frac{1}{1-\beta}} ]其中 (N) 是历史样本数(\beta) 是置信水平(C) 是与支撑集直径相关的常数。但在实际算例里我建议直接做参数扫描取 (\theta) 从0开始按步长0.01或者0.05逐步增大观察目标函数值总运行成本的变化趋势选择成本增长的“拐点”附近作为最终半径。这个拐点意味着再增大模糊集半径成本会大幅上升而鲁棒性提升有限。另外还有一个实用技巧如果数据集较小比如只有几十个历史样本半径的选择对结果非常敏感。这时我会采用留一法交叉验证即每次去掉一个历史样本用其余样本构造经验分布测试模型样本外表现选择平均表现最好的 (\theta)。2.4 Matlab中生成离散场景的完整流程我在Matlab里的场景生成流程大致如下读取历史误差数据err_hist维度是 (N \times d)(N) 为样本数(d) 为不确定参数个数风电、电负荷、热负荷等用kmeans或者同步回代消除法Backward Reduction将 (N) 个历史样本缩减为 (M) 个代表性场景计算每个场景的概率 (\pi_k n_k / N)其中 (n_k) 是第 (k) 个簇包含的原始样本数将缩减后的场景作为不确定参数的离散取值并在模糊集约束中引入概率不确定性。贴一段示意性Matlab代码% 假设 err_hist 为 N x d 的历史误差矩阵 rng(42); M 20; % 离散场景数 [idx, centroid] kmeans(err_hist, M, Replicates, 10); % 计算每个场景的经验概率 probs accumarray(idx, 1) / size(err_hist, 1); % 提取场景值每个聚类中心作为一个离散场景 scenarios centroid; % M x d % Wasserstein半径 theta 0.05; % 后续建模中每个场景k的误差取 scenarios(k,:)对应概率 probs(k)这里有几个细节需要注意kmeans本身是随机算法要设置Replicates避免陷入局部最优场景数 (M) 不宜太小否则会丢失历史数据中的尾部风险也不宜太大否则模糊集对偶后的变量数量成倍增长。我测试下来(M20\sim50) 对多数系统都适用同步回代消除法在Matlab里没有现成函数需要自己实现但效果往往比kmeans更适合风电场景因为它能保留极值场景。如果不会写可以用kmeans加后处理把原始样本中误差最大的几个点单独作为场景再对剩余样本聚类。3. 电热综合能源系统模型拆解与Matlab编码3.1 核心设备建模CHP、电锅炉、储热电热综合能源系统的特殊之处在于“耦合”。我最常用的设备建模方法如下。热电联产机组CHPCHP是电热耦合的核心。它同时产生电功率 (P_{chp}) 和热功率 (H_{chp})但两者不是独立的存在一个“热电比”可行域。简化的可行域可以写成[ P_{chp}^{\min} \leq P_{chp} \leq P_{chp}^{\max} ] [ 0 \leq H_{chp} \leq H_{chp}^{\max} ] [ H_{chp} \alpha P_{chp} \beta ](\alpha) 是热电比(\beta) 是常数项。实际中CHP的可行域更接近多边形可以用多组线性约束包围。在Matlab中用Yalmip定义时我会把每个线性约束单独写到cell数组中方便后续追加场景相关约束。电锅炉电锅炉本质是一个“电转热”设备消耗电功率 (P_{eb})产生热功率 (H_{eb})转换效率为 (\eta_{eb})[ H_{eb} \eta_{eb} P_{eb} ]它虽然没有复杂的耦合可行域但要注意启停次数约束和爬坡约束。在高比例风电场景下电锅炉常常被用来消纳弃风所以它的出力上限和爬坡速率会直接影响弃风量。储热罐储热罐是热力系统中的“电池”用储热量 (S_t) 表示[ S_{t1} S_t (H_{ch, t} - H_{dis, t}) \Delta t ]需要限制储热容量上下限、充放热功率上限以及同一个时刻不能同时充放热。这个“同时性”约束可以用二进制变量处理% 假设 H_ch 和 H_dis 是 sdpvar 变量 % b 是二进制变量1表示放热0表示充热 H_ch H_ch_max * b; H_dis H_dis_max * (1 - b);如果模型规模较大用这种大M法会引入二进制变量求解速度受影响可以先用线性互补约束或者把“同时充放”作为软约束看结果是否合理再决定是否加二进制变量。3.2 电力网络与热力网络的简化处理电力网络电力网络我用直流潮流模型节点功率平衡[ P_i \sum_{j \in \Omega_i} P_{ij} L_i ](P_{ij}) 是支路潮流与相角差成正比。在分布鲁棒优化中如果每个场景都要写一套全网络潮流约束变量数会爆炸。一个常用技巧是把网络约束预计算成机组出力向量到节点功率注入的映射矩阵也就是节点转移因子矩阵PTDF然后用矩阵乘法生成约束。在Matlab中可以用Matpower获取节点导纳矩阵再算PTDF。不过如果研究重点不在电力系统潮流也可以做更激进地简化忽略网络拓扑只保留节点功率平衡和线路备用容量约束。这样优化模型会小很多适合先跑通算法流程再逐步加网络约束。热力网络热力网络是最容易让人陷进去的地方。完整的热力管网模型包含供回水温度、管道流量、节点温度混合等非线性方程其中有一类双线性项流量×温度会让模型变成非凸的。如果多离散场景分布鲁棒再加上热网非线性模型求解难度会直线上升Gurobi直接罢工。我的经验是在电热综合能源系统优化中如果目标是“调度策略”而不是“管网水力热力详细仿真”完全可以用“节点热功率平衡储热/热惯性简化模型”来近似。也就是把热负荷看作节点功率需求热源通过热功率平衡方程供给忽略温度动态过程改为在约束中叠加一个“热惯性系数”允许短时热功率缺额在一定范围内用储热罐的热容量来吸收波动。这样虽然牺牲了热网温度分布的细节但整个模型保持线性分布鲁棒算法能够稳定收敛。如果你的课题确实需要精细热网模型建议把热网子问题单独拆出来用顺序迭代法先优化电/热功率分配再校验热网温度偏差过大则修正约束而不是把非线性直接塞进主优化模型。3.3 目标函数与决策变量目标函数一般是系统总运行成本最小包括CHP燃料成本与电出力和热出力都相关可用二次函数或分段线性近似向上级电网购电成本弃风惩罚成本储热充放损耗分布鲁棒优化中还会包含最坏分布下第二阶段调整成本。写成数学形式[ \min_{x} ; c(x) \sup_{P \in \mathcal{D}} \mathbb{E}_P \left[ Q(x,\xi) \right] ]其中 (x) 是第一阶段的“here-and-now”决策如机组启停、购电计划、储热充放计划(\xi) 是不确定参数风电出力、负荷误差(Q(x,\xi)) 是第二阶段根据实际场景做出的调整成本函数。在Matlab/Yalmip中我先定义第一阶段变量P_chp sdpvar(nt, n_chp); H_chp sdpvar(nt, n_chp); P_eb sdpvar(nt, 1); S sdpvar(nt1, 1); P_buy sdpvar(nt, 1); P_wind sdpvar(nt, 1); % 实际消纳风电再为每个离散场景定义第二阶段调整变量比如弃风量 (P_{curt}^{k})、切负荷量 (P_{shed}^{k})、设备出力调整量等。这些调整变量只在第二阶段场景中出现用于量化“预测值偏差带来的额外成本”。3.4 YalmipGurobi搭建优化模型的基本框架我用的是Yalmip做建模Gurobi做求解器。Yalmip的语法非常简洁适合快速搭建原型。下面是一个骨架% 定义变量第一阶段 P_chp sdpvar(nt, n_chp); ... % 定义场景变量第二阶段 P_curt sdpvar(nt, M); P_shed sdpvar(nt, M); % 约束 Constraints []; % 基础约束设备出力上下限、热/电功率平衡、储热动态 for t 1:nt Constraints [Constraints, P_chp(t,:) 0]; % ... 平衡约束 end % 场景相关约束第二阶段 for k 1:M for t 1:nt Constraints [Constraints, P_curt(t,k) 0]; % ... 风电实际出力与预测误差的关系 end end % 目标函数 Objective sum(CHP_cost) sum(P_buy_cost) ...; % 设置求解器 options sdpsettings(solver,gurobi,verbose,2); optimize(Constraints, Objective, options); % 后处理取第一阶段结果 P_chp_opt value(P_chp);注意Yalmip里不能把P_curt定义成三维变量时直接用于for循环的复杂操作尽量用矩阵运算。如果遇到Yalmip报“无法处理非凸约束”检查是不是出现了两个sdpvar相乘或绝对值/最大最小函数尽量用线性化改写。4. 多离散场景分布鲁棒优化的求解算法4.1 两阶段框架主问题与子问题如何衔接分布鲁棒优化的求解框架我采用经典的“两阶段鲁棒”思路第一阶段决策 (x) 必须在不确定性发生之前确定例如机组出力基点、购电合同量这些决策对未来所有场景都是一样的第二阶段决策 (y_k) 是看到场景 (k) 后做的调整例如弃风、切负荷、调整电锅炉功率用来保证功率平衡和约束不越限。如果直接用期望成本则模型为[ \min_{x} \left{ c^\top x \sum_{k1}^M \pi_k Q(x, \xi_k) \right} ]其中 (\xi_k) 是第 (k) 个离散场景(\pi_k) 是它的概率。这个形式看起来就像普通随机规划只要在每个场景分别写一遍 (Q(x, \xi_k)) 的约束把场景概率乘到目标函数里就可以整体求解。但这里是“分布鲁棒”也就是说场景概率 (\pi_k) 不是固定不变的而是可以在模糊集允许的范围内波动以最大化期望成本。所以真正的模型是[ \min_{x} \left{ c^\top x \max_{\pi \in \mathcal{D}} \sum_{k1}^M \pi_k Q(x, \xi_k) \right} ](\mathcal{D}) 是分布模糊集。好在对于Wasserstein模糊集内层的“max”可以通过线性规划对偶转换成一个“min”问题最终整体可以合并成一个大的MILP混合整数线性规划。这也是为什么我在前面强调一定要用Wasserstein模糊集因为它能让内层max问题有好的线性对偶形式。4.2 对偶转化Wasserstein模糊集如何变成线性约束这部分是算法能落地的核心我在Matlab里踩了很多次坑最终总结出一个相对稳健的转化方式。假设第二阶段问题 (Q(x, \xi)) 可以写成如下线性规划[ Q(x, \xi) \min_{y} ; q^\top y ] [ \text{s.t.} ; T y \geq h - M_x x - M_\xi \xi ]其中 (\xi) 是扰动参数。对于包含Wasserstein模糊集的分布鲁棒两阶段问题可以利用拉格朗日对偶把内层最坏期望转换为[ \min_{\lambda} ; \lambda \theta \frac{1}{N} \sum_{i1}^N \left( \text{某个与样本相关的子问题最优值} \right) ](\lambda) 是对偶变量(\theta) 是Wasserstein半径。子问题通常是一个线性规划其规模与历史样本数 (N) 和第二阶段变量数相关。听起来复杂但实际用Yalmip实现时我可以“手动”把这个对偶形式写成约束也可以采用更取巧的方法如果场景数量不是特别大直接用Yalmip的implies和sum写出max模型再让求解器内部分支定界处理。不过后者速度慢稍微大一点的系统就会卡住。我的建议是不要把对偶推导完全交给求解器而是自己先把对偶问题用公式化简好再进入Matlab编码。化简后的模型里每个历史样本对应一组辅助变量加入约束即可。典型的Yalmip代码结构如下% lambda为对偶变量 lambda sdpvar(1,1); % 每个历史样本对应一个子问题最优值变量 sub_val sdpvar(N,1); Constraints [Constraints, lambda 0]; for i 1:N % 根据子问题对偶互补条件写出 sub_val(i) 与变量 xi_i 的关系 Constraints [Constraints, sub_val(i) expr_with_xi_i]; end % 最坏期望成本近似为 lambda*theta mean(sub_val) DR_cost lambda * theta sum(sub_val) / N; Objective deterministic_cost DR_cost;这里最核心的是expr_with_xi_i的写法要依据你的系统模型来确定。如果模型里包含二进制变量第二阶段问题 (Q(x,\xi)) 不再是纯LP而是MILP那么对偶转化会非常麻烦通常需要引入二元变量表示互补条件。这种情况下我更推荐直接剖分场景把概率不确定放在一个辅助线性规划里用迭代算法例如列约束生成CCG来求解而不是一次性写出整个对偶模型。4.3 求解器选型和计算时间控制求解器方面首选Gurobi其次是Cplex最后才是Matlab自带的intlinprog。Gurobi对大规模MILP的求解速度非常快而且支持多线程在Yalmip中只需要设置options sdpsettings(solver,gurobi)即可。不过如果系统规模很大直接整体求解MILP仍然会面临“指数爆炸”。我实际测试过一个中等规模的算例10节点电网6节点热网20个离散场景第二阶段变量约1000个MILP有大约300个二进制变量Gurobi默认参数下需要十几分钟才能收敛到可行解这个速度对研究来说可以接受但对实时调度不够友好。加速办法有三个Benders或CCG分解把主问题和子问题分开求解主问题只包含第一阶段变量子问题逐场景求解并回传割平面。这种算法能显著减少同时求解的场景数但对编程能力要求高。场景裁剪不是每个历史样本都参与构造模糊集先用启发式方法剔除明显冗余的样本。比如聚类后只保留每个簇中心附近的若干样本减少对偶变量个数。松弛二进制变量如果启停变量不是关键先固定或松弛为连续变量求解一个LP作为初值再恢复二进制变量进行MIP搜索能大大减少分支次数。我的习惯是先跑一个“松弛LP”验证模型正确性再逐步恢复整数约束最后才考虑分解算法。这样能减少编程调试时间。4.4 线性化Max/Min/绝对值/双线性项怎么处理分布鲁棒优化中最容易出现的非线性项包括目标函数中的max、min、绝对值约束中的分段函数两个连续变量相乘如热网中的流量×温度二进制变量与连续变量相乘。这些在Yalmip中如果直接写虽然Yalmip能识别一部分但往往会把模型标记为“非凸”Gurobi可能不接。我的处理原则绝对值引入辅助变量 (u \geq |z|)等价于 (u \geq z) 且 (u \geq -z)max/min如上辅助变量加线性约束二进制×连续变量引入大M法变成两个线性约束例如 (0 \leq w \leq M b)(z - M(1-b) \leq w \leq z M(1-b))双线性项尽量避免实在避免不了就分段线性化或迭代求解。在电热综合能源系统中最容易忽略的是电功率平衡方程里风电出力与实际消纳量的关系。如果不允许弃风那就直接令风电出力等于预测值加误差如果允许弃风就要加一个“实际消纳风电”变量它不能超过可用风电也不能小于0。这个关系本身是线性的别多写乘积项。5. 算例结果与参数敏感性实战分析5.1 测试系统设置与参数为了验证算法我用了一个小型的电热综合能源系统配置如下参数数值电网节点数6热网节点数4CHP机组2台总电出力60MW热出力50MW电锅炉1台功率30MW效率0.9储热罐容量150MWh最大充放功率30MW风电装机50MW历史样本数500离散场景数20调度周期24小时负荷数据、风电预测数据来自公开数据集误差序列用历史残差拟合。基准场景为预测值误差场景通过kmeans聚类生成。模型中包含直流潮流约束和热功率平衡简化约束没有加完整热网水力模型。5.2 不同模糊集半径对调度成本的影响我做了 (\theta) 从0到0.5的扫描步长0.01结果如下(\theta)总运行成本万元弃风电量MWhCHP总出力MWh0随机规划18.622.49060.0519.219.89120.1020.116.59200.2022.311.29350.5025.85.3952可以明显看到随着 (\theta) 增大系统应对最坏分布的能力增强弃风量减少因为预留了更多消纳空间但总成本上升。在 (\theta 0.05\sim0.10) 之间成本上升幅度相对平缓超过0.1以后成本增长加速。我后来把 (\theta0.08) 作为基准值因为它能在不明显抬高运行成本的同时显著降低样本外最坏成本。这个结果也说明了一个道理模糊集半径并不是越大越好。半径过大意味着模型认为真实分布可能偏离历史经验分布非常远这相当于承认历史数据不可靠那还不如直接用纯鲁棒优化。5.3 离散场景数量的权衡我将离散场景数 (M) 从5增加到50设置 (\theta0.08)观察计算时间和目标函数值场景数 (M)目标函数值万元求解时间秒519.841019.6212019.2963019.13105019.01300目标函数值随着场景数增加而逐渐下降但边际收益越来越小求解时间却近似指数增长。20个场景已经能覆盖大部分分布信息50个场景成本只比20个低约1%求解时间却长了十几倍。因此我在论文和实际项目中都用20个场景作为默认值。这里要提醒一句场景数 (M) 和原始历史样本数 (N) 是两回事。(N) 决定模糊集中的样本基数(M) 是离散化后的支撑点数量。不要为减少计算量而把 (N) 也降到几十个否则 (N) 太小会导致经验分布本身不可靠(\theta) 再大也没用。5.4 常见报错与调试方法我在MatlabYalmipGurobi环境下遇到的问题基本可以归为几类第一类Yalmip报“无法验证凸性”通常是模型里出现了两个sdpvar相乘或者max/min/abs函数直接套在sdpvar上。先用check(Constraints)查看约束类型定位后手动线性化。第二类Gurobi报“Model is infeasible”不可行问题最常见的来源是第一阶段变量取值太紧导致第二阶段没有任何可行调整空间。我一般会在模型里加松弛变量比如允许微小切负荷和弃风并在目标函数加很大惩罚系数。先让问题可解再逐步减小松弛观察约束的紧度。也可以用ops sdpsettings(solver,gurobi,savesolveroutput,1)让Gurobi把不可行约束的 IIS不可行子系统导出来快速定位冲突约束。第三类求解器“Out of memory”这是场景数太多导致变量爆炸。解决办法减少离散场景数、去掉不必要的辅助变量、改用更紧凑的矩阵表示、把热网模型简化。很多情况下一个约束如果写成循环嵌套Yalmip会生成大量临时变量内存占用暴涨。尽量用矩阵运算一次定义整块约束。第四类结果出现奇怪的“零解”通常是目标函数系数写反了或者某个约束的下标索引从1开始但Matlab变量定义成了从0开始。这类问题最花时间没有捷径只能逐段检查。6. 我踩过的坑和给后来者的建议6.1 不要一上来就追求完整模型分步调试是王道我第一次做这个课题时想着一步到位把CHP、储热、电网、热网、分布鲁棒模糊集全部塞进模型结果代码写了三四百行报错信息铺天盖地。后来我把步骤拆成五级先做纯确定性优化也就是把风电预测值当作确定值验证设备模型和能量平衡约束正确加入多离散场景但暂时固定概率 (\pi_k)退化成随机规划验证场景生成和后处理逻辑正确引入模糊集先设 (\theta) 为很小的正数验证对偶转化后的模型不会爆炸逐步增大 (\theta)观察结果是否符合“越保守成本越高”的规律再精细化热网模型和网络拓扑。每一步都验证无误再进入下一步可以节省大量调试时间。6.2 热网模型不要过度精确双线性项是个大坑我在前面多次提到热网双线性项。这里再具体说一个案例我用标准节点法建模热网管道热功率 (Q m c_p (T_s - T_r))其中质量流量 (m) 和温差 ((T_s-T_r)) 都是变量乘积就是双线性。放到分布鲁棒优化中Gurobi报“Quadratic constraints not supported”还是好的有时候会直接报“Model is non-convex”。即使支持求解时间也会指数级上升。后来我换成了“质调节”模式假设供回水温度恒定管网中质量流量根据热负荷比例分配这样热网水力约束变成线性温度调节则由热源侧统一控制。这种简化在工程上完全说得通热网有热惯性短时间温度变化不大。研究重点应该放在多场景调度决策上而不是管网水力学细节。6.3 数据清洗比算法更重要异常值会毁掉模糊集有一段时间我只看模型不查数据结果用了一组包含明显坏数据的历史误差序列模糊集半径自动选到了0.4调度成本比纯鲁棒还要高。后来发现是历史数据里混了几天设备检修的记录误差值夸张到离谱。从那以后我每次构造模糊集之前都会先做异常值检测和分布可视化画误差直方图观察是否出现长尾用箱线图检查是否存在超过3倍四分位距的样本剔除异常值后重新聚类观察场景中心是否发生明显变化。这一步虽然简单但对最终结果的影响非常大。6.4 模糊集半径的选取不能只靠理论公式理论公式给出的是一个参考范围不是标准答案。我最初严格按照公式选 (\theta)结果模型过于保守。后来换成“成本-鲁棒性权衡曲线”来选反而更直观如果曲线斜率变化不大说明该点附近成本对鲁棒性提升的“性价比”较高。你甚至可以根据实际工程需求设定一个可接受的成本上限然后反推最大半径。6.5 后续可以怎么扩展这个算法框架其实有很强的扩展性。如果你想继续深入可以考虑把碳交易成本纳入目标函数因为电热综合能源系统减碳潜力很大碳价的不确定性可以用另一个模糊集建模引入需求响应尤其是热负荷的柔性可调能力让第二阶段的调整变量从“切负荷”变成“温度微调”更贴近实际把电转氢、电转气设备加入系统研究多能流耦合下的分布鲁棒调度将单时段优化扩展为多时段滚动优化用MPC的方式滚动更新模糊集能更好适应负荷波动。我在实际项目中已经尝试过把碳交易和需求响应加进去模型虽然更大但由于框架保持线性Gurobi依然能处理。如果遇到规模爆炸再考虑CCG分解不要一开始就上重型算法。最后说一点个人体会分布鲁棒优化看起来数学门槛高但真正落地到Matlab代码时核心工作其实是两件事——把不确定性建模成离散场景把对偶转化后的线性约束写对。前者靠数据清洗和场景缩减后者靠对线性规划的敏感和Yalmip的熟练使用。只要这两步踏踏实实做了后面的求解其实就是在Gurobi里“等结果”的过程。希望这篇文章能帮你在电热综合能源系统优化的路上少踩几个坑尤其是别像我一样在热网模型和异常数据上白白浪费两个星期。
