企业官网建设流程全解析

1. 先把名头拆开看：高热点算法背后的三个关键词

“高热点算法”这个说法第一次出现的时候，我以为是某个公众号的营销标题，但仔细想想，这个课题方向确实把近五年能源系统优化里最“卷”的三个关键词全凑齐了——数据驱动、分布鲁棒、电热综合能源系统。我第一次认真研究这套组合，是因为课题组要做一个含热电联产的区域综合能源系统调度模型，当时我还在用传统的随机优化，结果连续被审稿人追问“风电出力的分布函数是哪里来的”“假设正态分布的依据是什么”，这才逼着我把分布鲁棒优化这套方法完整走了一遍。

这个方向要解决的核心问题很直白：新能源出力预测不准，传统随机优化又太依赖我们人为假设的概率分布，一旦假设错了，最优解就跟着错。而纯鲁棒优化只盯着最坏情况，做出来的调度方案往往保守到没法用。分布鲁棒优化（DRO）走了一条中间路线——它不假设一个精确的概率分布，而是用历史数据构造一个“分布集合”，在这个集合里找最坏的那个分布来做决策。这套思路配合电热综合能源系统，正好能发挥电、热两种能源载体互相支撑的优势，所以最近几年论文发表量特别大，工程应用也在快速跟进。

这篇文章不打算堆公式，我会从工程实现的角度，把整套链路从头到尾讲透：电热IES设备怎么建模、多离散场景怎么用历史数据生成、模糊集怎么构造、min-max-min模型怎么一步步转成可求解形式，最后给出能跑的Matlab代码框架。适合正在做综合能源优化方向的学生、刚接触DRO的工程师，以及想把不确定性建模从随机优化升级到分布鲁棒的课题组成员参考。

1.1 随机优化和鲁棒优化为什么不够用

这三代方法的演进，本质上是一个“对不确定性信息掌握程度”的问题。随机优化假定我们完全知道风电出力、负荷等随机变量的精确概率分布，然后把期望成本降到最低。它的计算框架成熟，理论结果漂亮，但工程上有一个致命问题：真实分布我们根本不知道，通常只能假设一个正态分布或者用历史数据拟合。拟合出来的分布和真实分布之间的偏差，会让最优解在真实场景下表现很差，这种现象叫“分布误差导致的后悔值”。

鲁棒优化反过来走极端，它只要求知道随机变量的支撑集，比如“风电出力一定是0到80MW之间”，然后在这个区间内寻找最坏情况下的最优决策。它不依赖任何分布假设，结果非常稳妥，代价是过度保守——它把所有发生概率极小的极端事件和最可能发生的场景等权重对待，最终调度成本往往比随机优化高出10%甚至更多。对实际运行来说，这种方案意味着大量的弃风弃光、机组频繁调节，经济性很难接受。

分布鲁棒优化恰好补上了两者的空当。它的思想是：我不完全相信历史数据，但也承认历史数据能提供不少信息。于是用一个以经验分布为中心的“模糊集”来描述真实分布的取值范围，然后在这个模糊集内寻找使期望成本最大的“最坏分布”，并在此基础上做最优决策。这样既消化了数据信息，又保留了鲁棒性，代价只是多解一个凸优化问题。

1.2 分布鲁棒优化的核心思想：模糊集

为什么叫“模糊集”？它不是模糊数学那个概念，而是指“一系列候选概率分布的集合”。真实分布不知道，但我们相信它不会离经验分布太远，于是用一个几何球把这些可能的分布框住，球心就是历史数据得到的经验分布，球的半径控制保守程度。半径越大，集合越保守，决策越趋向于鲁棒优化；半径太小，集合接近单点，问题退化成随机优化。

构造模糊集的方式有好几种。最早的是矩模糊集，约束候选分布的均值、协方差落在某个范围内，理论推导起来方便，但实际取矩约束的边界参数很困难。这几年最火的是基于Wasserstein距离的模糊集，它跟概率分布之间的距离有关，有比较好的理论保证和有限样本界的结论，对偶转化也比较漂亮。还有一类在中文文献里特别常见的做法，是用1-范数和∞-范数去约束场景概率向量和名义概率之间的偏差，因为转化之后是线性约束，工程实现最简单，正好跟“多离散场景”天然匹配。本文代码实现选择的就是这种范数约束型模糊集。

1.3 电热综合能源系统为什么特别适合这套框架

电热综合能源系统是一个物理结构复杂、不确定性传导路径多的系统。风电、光伏的波动直接冲击电力平衡，而热负荷波动、管网延迟、热惯性又从另一个时间尺度影响系统运行。CHP机组把电和热耦合在一起，电出力的调节必须跟着热出力走，风电上来了电不好消纳，就只能通过电锅炉把多余的电转成热，蓄热罐再把热能存起来。这一套调节链条里，每一步都涉及不确定性，不是单纯的“电力平衡”问题。

如果只做确定性优化，风电预测误差会被后备机组吃掉，调度方案谈不上最优。如果用分布鲁棒框架，我们能把历史运行数据直接用于模糊集构造，让系统在最坏但合理的分布场景下，寻找成本和风险之间的平衡点。再加上热网本身具备天然的储能特性和时间延迟，能在阶段间提供额外的灵活性，使得两阶段分布鲁棒模型的第二段决策有更大调整空间。这也是为什么最近几年关于“数据驱动分布鲁棒+电热IES”的论文会大量出现的原因——问题本身、数据条件、方法工具三者确实是匹配的。

2. 电热综合能源系统建模：从设备机理到优化约束

模型是整个优化问题的基础。我直接把常用的区域级电热IES模型拆开讲，这个模型覆盖了最典型的设备：抽凝式CHP机组、电锅炉、蓄热罐、风电场、普通电负荷和热负荷，外加一个简化的直流潮流电网和一个热力节点平衡热网。设备级模型不搞复杂机理，全部用线性约束描述，因为只有线性模型才能在CCG框架里高效求解。

2.1 CHP机组可行域：电热耦合的核心

CHP机组是整个系统“电热耦合”最关键的一台设备。抽凝式CHP的电出力和热出力不是完全独立的，它们共同落在某个二维可行域内。常见做法把这个可行域近似成凸多边形，用一组线性不等式来刻画。典型约束如下：

$$ P_{i}^{\min} \le P_{i,t}^{e} \le P_{i}^{\max} $$

$$ H_{i}^{\min} \le H_{i,t}^{h} \le H_{i}^{\max} $$

$$ P_{i,t}^{e} + \lambda_{i} H_{i,t}^{h} \le \Gamma_{i} $$

第一条和第二条是电、热出力的上下限，第三条是电热耦合限制，$\lambda_i$ 表示抽汽系数，$\Gamma_i$ 反映汽轮机进汽量的上限。实际建模型时，视精度要求还可以再加几条不同斜率的约束把可行域刻得更细，有些文献直接把CHP运行区间画成一个梯形、四边形甚至五边形，本质都是这一条条线性不等式拼出来的。

值得提醒的是，背压式CHP是抽凝式的一个极端特例，电出力和热出力严格成正比，相当于可行域退化成一条线段，建模更简单但灵活性差很多。目前绝大多数IES调度论文用的是抽凝式，因为留了调节空间，问题更有挑战性也更有研究价值。

2.2 电锅炉、蓄热罐与热网平衡

电锅炉是电能转热能的装置，模型非常简单：

$$ H_{t}^{EB} = \eta_{EB} \cdot P_{t}^{EB} $$

$P_t^{EB}$ 是消耗的电功率，$H_t^{EB}$ 是产出的热功率，$\eta_{EB}$ 是电热转换效率，一般取0.9到0.98之间。电锅炉是风电消纳的利器——当风电大发、电网平衡出现困难时，把多余的电转成热供给热负荷或蓄热罐，能有效缓解弃风。

蓄热罐的模型核心是能量状态递推：

$$ SOC_{t+1} = SOC_t + \eta_{ch} P_{t}^{ch} - \frac{P_{t}^{dis}}{\eta_{dis}} $$

$SOC_t$ 是蓄热量，$P_t^{ch}$ 和 $P_t^{dis}$ 分别是充放热功率，$\eta_{ch}$ 和 $\eta_{dis}$ 是充放热效率。除此之外还需要约束充放热功率上下限、蓄热容量上下限，以及严禁同时充放热。蓄热罐的作用可以类比成一个热水“充电宝”，它把热力系统的时间灵活性释放出来，让CHP机组不必在某个时刻硬顶着热负荷运行，而可以根据电价或风电情况灵活调节电出力。

热网层面，如果不想陷入管道水力瞬态建模的泥潭，区域级调度模型一般直接用热力节点平衡加简化网络约束：

$$ \sum_{source} H_{src,t} + \sum_{storage,dis} H_{dis,t} = \sum_{load} H_{load,t} + \sum_{storage,ch} H_{ch,t} $$

即所有热源产生的热量加上蓄热罐的放热量，等于热负荷消耗量加上蓄热罐的充热量。如果要做更精细的分析，需要引入供水/回水温度、管道传输延迟和热损耗模型，那一套属于详细热网动态仿真，和优化调度的尺度不完全一样。我的建议是：跑优化先用节点平衡模型，等找到调度规律之后再单独做热网水力校核。

2.3 直流潮流与电热耦合约束的衔接

电力网络部分我采用直流潮流模型，这是IES优化里的标准简化手段。直流潮流假设电压幅值近似恒定、支路两端相角差很小，只考虑有功功率：

$$ P_{k,t} = \sum_{j \in N} B_{kj}(\theta_{k,t} - \theta_{j,t}) $$

其中 $B_{kj}$ 是节点电纳矩阵的负元素，$\theta$ 是相角。节点功率平衡则写成：

$$ \sum_{g \in G_k} P_{g,t} + P_{k,t}^{wind} + P_{k,t}^{buy} = \sum_{l \in L_k} P_{l,t} + P_{k,t}^{EB} + \sum_{j} B_{kj}(\theta_{k,t} - \theta_{j,t}) $$

这条式子的物理意义是：注入节点的功率等于流出节点的功率。其中风电、购电是注入项，电负荷、电锅炉消耗是流出项。注意电锅炉在这里以负荷形式出现在等号右侧，它的耗电量和热力侧产热量通过效率公式连在一起，于是电网络和热网络就通过CHP和电锅炉这两个“交叉设备”完成了耦合。

电网的线路潮流约束也别忘了加：

$$

• P_{kj}^{\max} \le B_{kj}(\theta_{k,t} - \theta_{j,t}) \le P_{kj}^{\max} $$

这套约束在Matlab里其实就是一组标准的线性不等式，YALMIP或者直接写矩阵都可以处理，关键是建立节点、线路、发电机的索引关系，不要出现节点类型对不上或者支路电纳符号错误的问题。

2.4 目标函数和不确定性参数的定位

目标函数直接决定调度的“价值取向”。我常做的最小化总成本目标函数包含五个部分：

$$ \min \sum_{t} \left( C_{t}^{buy} P_{t}^{buy} + \sum_{i} f_i(P_{i,t}^{e}, H_{i,t}^{h}) + \lambda_{wind} \Delta P_{t}^{wind} + \lambda_{load} \Delta P_{t}^{load} \right) $$

第一项是向外部电网买电的成本，$C_t^{buy}$ 是分时电价；第二项是CHP机组的燃料成本，一般写成关于电出力和热出力的二次函数，在优化里常用分段线性近似；第三项是弃风惩罚，$\Delta P_t^{wind}$ 是被迫削减的风电功率，惩罚系数 $\lambda_{wind}$ 设得很高，表达“尽量消纳风电”的意愿；第四项是切负荷惩罚，正常情况下这个值应该严格为0，一旦出现说明系统在某个场景下确实存在失负荷风险，这也是分布鲁棒优化重点关注的风险来源。

不确定性参数如何定位也很关键。最常见的选择是把风电出力作为随机变量，因为它的预测误差最大，对电力平衡的冲击最明显。更进一步可以同时把电负荷、热负荷的预测误差也划入随机向量，但注意随机变量维度升高之后，场景生成和模糊集参数标定的难度都会增大，初学阶段建议先只做风电随机，模型跑通了再扩展维度。我在后面的代码实现里，就默认只把风电出力视为随机变量。

3. 多离散场景怎么来：历史数据驱动的场景生成与缩减

标题里“数据驱动”四个字不是噱头，它在实现上就体现在两个地方：一是用历史数据生成多离散场景，二是用历史数据构造模糊集。这一章先解决场景的来源问题。场景生成的目标是把连续型的随机变量转化为有限个离散场景，每个场景带一个概率，最终形成一个S个场景的离散概率分布。

3.1 随机变量辨识与数据处理

第一步不是直接生成场景，而是先搞清楚历史数据里到底哪些量能反映不确定性。对风电场，典型的做法是收集一段时间的实际出力数据和对应的预测数据，然后把误差序列提取出来：$\xi_t = P_{t}^{actual} - P_{t}^{forecast}$。这个误差序列就是建模随机变量的原始素材。为什么不直接建模风电出力？因为出力本身和时间、季节强相关，规律太复杂，而预测误差更接近零均值、相对平稳，对分布假设的鲁棒性更强。

数据量方面，至少要有几个月的数据才谈得上“数据驱动”。如果只有几天的数据，构造出来的经验分布和模糊集都缺乏统计意义，那还不如直接用鲁棒优化。另外，风电功率有天然的物理边界，场景生成之后一定要做边界检查，把所有小于0或大于装机容量的数值clip到合理区间，否则优化模型会给出完全违背物理的调度结果。

3.2 蒙特卡洛批量生成场景

拿到历史误差序列后，常见的做法是用蒙特卡洛抽样生成大量初始场景。最简单的抽样方式是Bootstrap重采样：从历史误差样本里随机抽取一天或多个时段，组合成一个完整的24小时误差场景，然后叠加到当天的预测曲线上。这种方法的好处是保持了误差之间的时序相关性，不会出现那种“前一个小时误差+30%、后一个小时误差-30%”的过度振荡场景。如果数据量充分，也可以用核密度估计拟合误差的分布函数，再用函数抽样，但那样引入的分布假设又会掩盖一部分“数据驱动”的初衷。

Matlab代码可以这样写：

% wind_hist: 历史风电出力矩阵，大小 N x M_hist（N为时段数，M_hist为历史天数） % wind_forecast: 当前预测风电出力，N x 1 % M_total: 需要生成的初始场景数 err_hist = wind_hist - wind_forecast; % 近似误差样本 wind_scenes = zeros(N, M_total); for s = 1:M_total idx = randi(size(err_hist,2), 1); wind_scenes(:, s) = wind_forecast + err_hist(:, idx); end wind_scenes(wind_scenes < 0) = 0; wind_scenes(wind_scenes > wind_cap) = wind_cap;

代码里需要注意，err_hist(:, idx)取出的是同一天的整段误差序列，这比逐时段独立抽取更接近真实时序特性。生成500到1000个初始场景后，下一步就是缩减。

3.3 聚类缩减：谁留下，谁合并

初始场景可能有上千个，直接塞进两阶段优化模型里，求解时间会增长到难以接受。实际工程里一般用聚类算法把场景数缩减到10到20个，再把这些场景作为离散代表场景参与优化。K-means是最直接的选择，原因在于接口成熟、速度极快、效果稳定。

% 把 N x M_total 的矩阵转置成 M_total x N 再聚类 [idx, C] = kmeans(wind_scenes', S, 'Start', 'plus', 'MaxIter', 500); % 计算每个聚类的场景概率 prob_s = histcounts(idx, S)' / M_total; % 簇中心作为代表场景 wind_rep = C'; % N x S 的代表场景矩阵

有一个工程细节要特别说明：K-means聚出来的簇中心是各簇的质心，它未必落在原始样本点上，有时会出现轻微越界或破坏边界条件。我在代码里会对簇中心再次做物理边界clip。另外，如果历史样本分布非常不均匀，K-means可能会产生个别样本极少的簇，代表场景的概率接近0，这种情况要手动合并到最近的大簇里，否则那个场景对应的分布概率在模糊集里会出现数值病态。

3.4 场景数量与概率分布的工程选择

场景数量S不是越大越好，也不是越小越好，它直接影响两阶段模型的求解规模和模糊集的保守程度。以我的经验，初始生成500个场景、聚类缩减到15个左右，是一个计算时间和精度都合理的平衡点。S=10以下，分布信息丢失明显，DRO的模糊集再大也无法弥补信息不足；S超过30，CCG主问题的规模会膨胀，求解时间从几分钟涨到几十分钟，收益却微乎其微。

场景概率的初始分布 $p_s^0$ 我直接用聚类簇的频率占比：$p_s^0 = n_s / M_{total}$，其中 $n_s$ 是第s个簇包含的样本数。这样得到的 $p_s^0$ 也叫做名义分布，它是模糊集的“球心”。后面构造模糊集时，所有允许的候选分布都围绕这个名义分布波动。

4. 分布鲁棒模型构造与求解转化：从min-max-min到可计算形式

场景有了，模型有了，接下来是最核心也最容易卡住的一步：把两阶段分布鲁棒问题写成数学上可求解的形式。很多初学者卡在这一步，并不是因为推导太难，而是没有把“max-min”嵌套结构里每一层到底在做什么想清楚。

4.1 基于范数约束的模糊集构造

我用的是范数约束型模糊集，原因在前面说过：它简单、线性、和离散场景天然契合。模糊集写成这样：

$$ \Psi = \left{ \mathbf{p} \in \mathbb{R}^{S} : \sum_{s=1}^{S} |p_s - p_s^{0}| \le \theta_1, \quad \max_{s} |p_s - p_s^{0}| \le \theta_{\infty}, \quad \sum_{s=1}^{S} p_s = 1, \quad p_s \ge 0 \right} $$

第一个约束是1-范数约束，控制所有场景概率相对名义分布的整体偏差，$\theta_1$ 越大，分布整体偏移的程度越强；第二个约束是无穷范数约束，限制单个场景概率的最大偏移，防止某个场景概率被过度抬高。这两个参数搭配使用，既允许整体分布“漂移”，又不让某一个场景单独变成主导。

参数 $\theta_1$ 和 $\theta_{\infty}$ 怎么取？理论上可以根据历史样本数量和置信度推导，工程上则建议通过灵敏度分析来标定。我的经验是：S=15时，$\theta_1$ 从0.3取到1.0之间，$\theta_\infty$ 从0.05取到0.2之间，一般来说 $\theta_\infty$ 是 $\theta_1$ 的1/5到1/10。$\theta_1=0$ 时模糊集退化为单点，模型退化成随机优化；$\theta_1$ 越来越大时，结果越来越保守，逐渐逼近纯鲁棒优化的效果。

4.2 两阶段分布鲁棒模型的标准形式

两阶段分布鲁棒问题可以写成标准的min-max-min结构：

$$ \min_{x \in X} ; c^T x + \max_{p \in \Psi} \sum_{s=1}^{S} p_s Q(x, \xi_s) $$

第一层的 $\min$ 是第一阶段决策，也就是在不确定性未实现之前就要确定的计划，比如外购电计划、蓄热罐的基准供热量。第二层的 $\max$ 是在模糊集内寻找最坏分布，它决定“如果真实分布是模糊集里最不利的那个，我的期望成本会是多少”。最内层的 $Q(x, \xi_s)$ 是第二阶段优化，给定第一阶段变量和某个具体的随机场景，决策系统内各可调设备的出力，最小化该场景下的运行成本。

第二阶段问题本质是一个线性规划：

$$ Q(x, \xi_s) = \min_{y \in Y(x, \xi_s)} ; d^T y $$

其中 $y$ 是第二阶段变量，包括各场景下CHP的出力调节量、电锅炉功率、蓄热罐充放热功率等；$Y(x, \xi_s)$ 是第二阶段可行域，依赖第一阶段变量和场景参数。整个结构的直观理解是：先定一个“计划”，再假设老天爷用模糊集里最坏的概率分布产生风电偏差，然后在这个最坏分布下用最优调度去兜底。

4.3 第二阶段的对偶处理与子问题变形

max-min嵌套求解的核心手段是把内层min问题转成对偶max问题，从而把“max-min”变成“max-max”，二者合并后依然是一个凸问题。具体来说，对每个场景的内层LP，利用强对偶条件写出它的对偶问题：

$$ Q(x, \xi_s) = \max_{\lambda \in \Lambda_s} ; \lambda^T (b_s - A_s x) $$

外层的 $\max_p \sum p_s Q(x, \xi_s)$ 结合每个 $Q(x,\xi_s)$ 的对偶形式，最终会变成一个带耦合约束的max问题。虽然推导过程看着复杂，但从工程实现角度，其实不用手动做所有对偶推导——只要保证内层子问题满足强对偶条件（原问题可行且有界），可以直接调用求解器求解，然后用KKT条件提取对偶乘子或直接计算割平面。

我在实际代码里采用了一种更工程化的做法：在CCG迭代的每一步，给定当前第一阶段解 $x^$，先对每个场景单独求解第二阶段LP，得到 $Q_s(x^)$；然后在外层求解一个只有S+2个变量的线性规划——最坏分布问题：

$$ \max_{p \in \Psi} ; \sum_{s} p_s Q_s(x^*) $$

这个外层LP极小，只有S个变量和S+4条约束，求解是微秒级的事。它得到的最坏分布 $p^*$ 直接将用于生成主问题的割平面。

4.4 为什么选CCG而不是Benders分解

解决两阶段分布鲁棒问题的经典算法有Benders分解和列与约束生成（CCG）两种。很多教材先讲Benders，导致初学者默认Benders是首选，但从实际收敛速度看，CCG在这类问题上明显更优。Benders通过对偶乘子添加割平面去逼近第二阶段价值函数，本质上是用一阶近似逐步逼近，需要多次迭代才能收敛；而CCG每一次迭代把最坏分布对应的场景决策变量直接复制进主问题，相当于把“阶梯”直接搭到真值附近，收敛速度要快得多。

CCG的直观理解是：每一次迭代都把最坏的那个分布“定格”下来，然后要求主问题在这个固定分布下优化。由于主问题拿到的割是原问题的下界，而子问题给出的期望成本是上界，上下界不断逼近，直到gap小于容差。对一个场景数S=15的中等规模IES模型，CCG通常在3到8轮迭代内收敛，而Benders往往要几十轮。这就是我选择CCG的根本原因。

5. Matlab代码实现：核心模块与关键代码解读

理论部分讲完了，这一章直接上代码。我实现的这套代码用YALMIP作为建模语言，求解器可以选Gurobi或者Cplex，两者的Matlab接口都成熟稳定。整体代码结构按照“数据准备—场景生成—主问题—子问题—CCG循环—结果后处理”的流程组织。

5.1 整体代码架构

整个项目我习惯分成以下几个文件，清晰易维护：

main.m % 主程序，控制CCG迭代流程 case_data.m % 系统参数定义：负荷、风电、设备参数、分时电价 gen_scenarios.m % 历史数据生成初始场景 + K-means聚类缩减 build_ambiguity.m % 构造模糊集参数 build_MP.m % 构建主问题（YALMIP模型） solve_SP.m % 子问题求解：内层场景LP + 最坏分布LP plot_result.m % 结果可视化

main.m的流程就三件事：加载系统参数、生成场景、启动CCG循环。参数全部集中在case_data.m里，换算例时只改这个文件即可。场景生成和模糊集构造各自独立成函数，方便单独调试。

5.2 场景生成模块的实现

场景生成模块我在第3章已经给出了核心代码片段，这里补充一点细节。在实际代码里，gen_scenarios.m函数返回三样东西：代表场景矩阵wind_rep（大小 N x S）、名义概率向量p0（长度 S）、以及初始聚类索引idx。其中idx在后续判断“哪些场景物理上相近”时会用到。

function [wind_rep, p0, prob_cluster] = gen_scenarios(wind_forecast, wind_hist, M_total, S, wind_cap) N = length(wind_forecast); err_hist = wind_hist - wind_forecast;