电热综合能源系统分布鲁棒优化：数据驱动模糊集与Matlab实现全解析

做电热综合能源系统优化有一段时间了，说实话，最头疼的不是模型本身，而是不确定性。天气一变，风电出力和热负荷预测误差就会被放大，这时候如果还沿用传统的确定性调度，一个极端场景就可能让系统直接越限。我最近把一套基于数据驱动+多离散场景分布鲁棒的电热综合能源系统优化算法，在Matlab里完整跑通了，从模糊集构造、场景生成到两阶段求解，踩了不少坑，也整理出不少可以复用的经验。这篇文章打算把整个思路从头到尾捋一遍：为什么传统的随机规划和鲁棒优化不够用、数据驱动模糊集怎么在Matlab里落地、电热综合能源系统模型怎么拆、两阶段求解框架怎么实现，最后附上算例调试和避坑心得。适合正在做电热综合能源优化、综合能源系统调度、分布鲁棒优化相关研究的同学参考，尤其适合那些已经有Matlab+Yalmip基础，但不知道如何把分布鲁棒算法落到具体算例里的朋友。

1. 为什么我放弃了传统的随机优化和鲁棒优化

1.1 电热综合能源系统的“不确定性”到底来自哪里

电热综合能源系统不是简单的电网加热网，而是通过热电联产机组、电锅炉、热泵、储热等设备把电力系统和热力系统耦合在一起。风电、光伏的间歇性让电力侧的源侧不确定性变得非常突出，而热负荷又和天气、用户行为强相关，尤其是在北方采暖季，昼夜温差会导致热负荷预测偏差明显。还有一类不确定性来自市场价格，如果园区需要从上级电网购电，现货电价波动也会直接影响调度决策。

这些不确定性不是“可忽略的小噪声”，而是可能达到预测值20%～30%的偏差。以风电为例，一台风机在调度时段的实际出力可能会低于预测出力一半以上。热负荷虽然没有那么剧烈，但管网的热惯性会让温度变化滞后，如果调度方案在热负荷预测基础上没有任何裕度，很容易出现末端用户供热不足。

我在最早做这个课题时，第一反应是用蒙特卡洛抽样生成大量场景，然后做期望值优化，也就是经典的随机规划。但实际跑下来发现两个问题：第一，随机规划非常依赖“概率分布已知”这个假设。实际工程里我们只有历史数据，不知道真实分布，一旦用历史经验分布代替真实分布，优化出来的方案看起来样本外表现不错，但遇到没见过的极端天气，照样翻车。第二，为了让期望成本有意义，通常需要生成上千个场景，模型规模急剧膨胀，Matlab里光构建约束矩阵就够吃内存，加上电热耦合约束，求解时间经常按小时算。

1.2 鲁棒优化“太保守”的根源在于集合构造

后来我又试了传统鲁棒优化，思路很简单：不确定参数在一个盒式集合或椭球集合里变化，我找到最坏情况下成本最小的方案。这种方法的优点是模型一般可以转化成确定性规划，求解快，而且理论上保证了任意的不确定实现下都可行。但代价也很明显——成本高得离谱。

举个例子，风电出力预测是100MW，盒式集合取±30%，鲁棒优化会认为风电可能只有70MW，然后让火电机组预留30MW的备用容量。实际上一百个历史样本里，可能只有一两个样本出力低于75MW，但鲁棒优化不为所动，坚持按最坏情况布置备用。这种保守性在纯电力系统里或许还能接受，放到电热综合能源系统里就会被放大，因为热力系统惯性大、允许短时间内有一定温度偏差，如果也让储热罐时刻保持“最坏情况下的最大容量”，整个系统的经济性会非常难看。

关键是，传统鲁棒优化完全不利用历史数据的分布信息。它甚至不关心“最坏情况发生的概率是多少”，只要落在集合内部，就按等概率处理。这在数据丰富的工程场景里明显是浪费。

1.3 分布鲁棒优化：在不确定分布中找“最坏”的期望

分布鲁棒优化（Distributionally Robust Optimization，DRO）的思路恰好落在两者之间。它假设不知道真实分布，但知道真实分布一定落在某一个由历史数据构造的模糊集（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里的场景生成流程大致如下：

1. 读取历史误差数据err_hist，维度是 (N \times d)，(N) 为样本数，(d) 为不确定参数个数（风电、电负荷、热负荷等）；
2. 用kmeans或者同步回代消除法（Backward Reduction）将 (N) 个历史样本缩减为 (M) 个代表性场景；
3. 计算每个场景的概率 (\pi_k = n_k / N)，其中 (n_k) 是第 (k) 个簇包含的原始样本数；
4. 将缩减后的场景作为不确定参数的离散取值，并在模糊集约束中引入概率不确定性。

贴一段示意性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) 不宜太小，否则会丢失历史数据中的尾部风险；也不宜太大，否则模糊集对偶后的变量数量成倍增长。我测试下来，(M=20\sim50) 对多数系统都适用；
• 同步回代消除法在Matlab里没有现成函数，需要自己实现，但效果往往比kmeans更适合风电场景，因为它能保留极值场景。如果不会写，可以用kmeans加后处理：把原始样本中误差最大的几个点单独作为场景，再对剩余样本聚类。

3. 电热综合能源系统模型拆解与Matlab编码

3.1 核心设备建模：CHP、电锅炉、储热

电热综合能源系统的特殊之处在于“耦合”。我最常用的设备建模方法如下。

热电联产机组（CHP）

CHP是电热耦合的核心。它同时产生电功率 (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_{t+1} = 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(nt+1, 1); P_buy = sdpvar(nt, 1); P_wind = sdpvar(nt, 1); % 实际消纳风电

再为每个离散场景定义第二阶段调整变量，比如弃风量 (P_{curt}^{k})、切负荷量 (P_{shed}^{k})、设备出力调整量等。这些调整变量只在第二阶段场景中出现，用于量化“预测值偏差带来的额外成本”。

3.4 Yalmip+Gurobi搭建优化模型的基本框架

我用的是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_{k=1}^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_{k=1}^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_{i=1}^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，那么对偶转化会非常麻烦，通常需要引入二元变量表示互补条件。这种情况下，我更推荐直接剖分场景，把概率不确定放在一个辅助线性规划里，用迭代算法（例如列约束生成C&CG）来求解，而不是一次性写出整个对偶模型。

4.3 求解器选型和计算时间控制

求解器方面，首选Gurobi，其次是Cplex，最后才是Matlab自带的intlinprog。Gurobi对大规模MILP的求解速度非常快，而且支持多线程，在Yalmip中只需要设置options = sdpsettings('solver','gurobi')即可。

不过，如果系统规模很大，直接整体求解MILP仍然会面临“指数爆炸”。我实际测试过一个中等规模的算例：10节点电网+6节点热网，20个离散场景，第二阶段变量约1000个，MILP有大约300个二进制变量，Gurobi默认参数下需要十几分钟才能收敛到可行解，这个速度对研究来说可以接受，但对实时调度不够友好。

加速办法有三个：

1. Benders或C&CG分解：把主问题和子问题分开求解，主问题只包含第一阶段变量，子问题逐场景求解并回传割平面。这种算法能显著减少同时求解的场景数，但对编程能力要求高。
2. 场景裁剪：不是每个历史样本都参与构造模糊集，先用启发式方法剔除明显冗余的样本。比如聚类后只保留每个簇中心附近的若干样本，减少对偶变量个数。
3. 松弛二进制变量：如果启停变量不是关键，先固定或松弛为连续变量，求解一个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 测试系统设置与参数

为了验证算法，我用了一个小型的电热综合能源系统，配置如下：

负荷数据、风电预测数据来自公开数据集，误差序列用历史残差拟合。基准场景为预测值，误差场景通过kmeans聚类生成。模型中包含直流潮流约束和热功率平衡简化约束，没有加完整热网水力模型。

5.2 不同模糊集半径对调度成本的影响

我做了 (\theta) 从0到0.5的扫描，步长0.01，结果如下：

可以明显看到，随着 (\theta) 增大，系统应对最坏分布的能力增强，弃风量减少（因为预留了更多消纳空间），但总成本上升。在 (\theta = 0.05\sim0.10) 之间，成本上升幅度相对平缓，超过0.1以后成本增长加速。我后来把 (\theta=0.08) 作为基准值，因为它能在不明显抬高运行成本的同时，显著降低样本外最坏成本。

这个结果也说明了一个道理：模糊集半径并不是越大越好。半径过大意味着模型认为真实分布可能偏离历史经验分布非常远，这相当于承认历史数据不可靠，那还不如直接用纯鲁棒优化。

5.3 离散场景数量的权衡

我将离散场景数 (M) 从5增加到50，设置 (\theta=0.08)，观察计算时间和目标函数值：

目标函数值随着场景数增加而逐渐下降，但边际收益越来越小；求解时间却近似指数增长。20个场景已经能覆盖大部分分布信息，50个场景成本只比20个低约1%，求解时间却长了十几倍。因此我在论文和实际项目中都用20个场景作为默认值。

这里要提醒一句：场景数 (M) 和原始历史样本数 (N) 是两回事。(N) 决定模糊集中的样本基数，(M) 是离散化后的支撑点数量。不要为减少计算量而把 (N) 也降到几十个，否则 (N) 太小会导致经验分布本身不可靠，(\theta) 再大也没用。

5.4 常见报错与调试方法

我在Matlab+Yalmip+Gurobi环境下遇到的问题，基本可以归为几类：

第一类：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、储热、电网、热网、分布鲁棒模糊集全部塞进模型，结果代码写了三四百行，报错信息铺天盖地。后来我把步骤拆成五级：

1. 先做纯确定性优化，也就是把风电预测值当作确定值，验证设备模型和能量平衡约束正确；
2. 加入多离散场景，但暂时固定概率 (\pi_k)，退化成随机规划，验证场景生成和后处理逻辑正确；
3. 引入模糊集，先设 (\theta) 为很小的正数，验证对偶转化后的模型不会爆炸；
4. 逐步增大 (\theta)，观察结果是否符合“越保守成本越高”的规律；
5. 再精细化热网模型和网络拓扑。

每一步都验证无误，再进入下一步，可以节省大量调试时间。

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依然能处理。如果遇到规模爆炸，再考虑C&CG分解，不要一开始就上重型算法。

最后说一点个人体会：分布鲁棒优化看起来数学门槛高，但真正落地到Matlab代码时，核心工作其实是两件事——把不确定性建模成离散场景，把对偶转化后的线性约束写对。前者靠数据清洗和场景缩减，后者靠对线性规划的敏感和Yalmip的熟练使用。只要这两步踏踏实实做了，后面的求解其实就是在Gurobi里“等结果”的过程。希望这篇文章能帮你在电热综合能源系统优化的路上少踩几个坑，尤其是别像我一样在热网模型和异常数据上白白浪费两个星期。