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_{t1} 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_histN为时段数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个左右是一个计算时间和精度都合理的平衡点。S10以下分布信息丢失明显DRO的模糊集再大也无法弥补信息不足S超过30CCG主问题的规模会膨胀求解时间从几分钟涨到几十分钟收益却微乎其微。场景概率的初始分布 $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_{s1}^{S} |p_s - p_s^{0}| \le \theta_1, \quad \max_{s} |p_s - p_s^{0}| \le \theta_{\infty}, \quad \sum_{s1}^{S} p_s 1, \quad p_s \ge 0 \right} $$第一个约束是1-范数约束控制所有场景概率相对名义分布的整体偏差$\theta_1$ 越大分布整体偏移的程度越强第二个约束是无穷范数约束限制单个场景概率的最大偏移防止某个场景概率被过度抬高。这两个参数搭配使用既允许整体分布“漂移”又不让某一个场景单独变成主导。参数 $\theta_1$ 和 $\theta_{\infty}$ 怎么取理论上可以根据历史样本数量和置信度推导工程上则建议通过灵敏度分析来标定。我的经验是S15时$\theta_1$ 从0.3取到1.0之间$\theta_\infty$ 从0.05取到0.2之间一般来说 $\theta_\infty$ 是 $\theta_1$ 的1/5到1/10。$\theta_10$ 时模糊集退化为单点模型退化成随机优化$\theta_1$ 越来越大时结果越来越保守逐渐逼近纯鲁棒优化的效果。4.2 两阶段分布鲁棒模型的标准形式两阶段分布鲁棒问题可以写成标准的min-max-min结构$$ \min_{x \in X} ; c^T x \max_{p \in \Psi} \sum_{s1}^{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^)$然后在外层求解一个只有S2个变量的线性规划——最坏分布问题$$ \max_{p \in \Psi} ; \sum_{s} p_s Q_s(x^*) $$这个外层LP极小只有S个变量和S4条约束求解是微秒级的事。它得到的最坏分布 $p^*$ 直接将用于生成主问题的割平面。4.4 为什么选CCG而不是Benders分解解决两阶段分布鲁棒问题的经典算法有Benders分解和列与约束生成CCG两种。很多教材先讲Benders导致初学者默认Benders是首选但从实际收敛速度看CCG在这类问题上明显更优。Benders通过对偶乘子添加割平面去逼近第二阶段价值函数本质上是用一阶近似逐步逼近需要多次迭代才能收敛而CCG每一次迭代把最坏分布对应的场景决策变量直接复制进主问题相当于把“阶梯”直接搭到真值附近收敛速度要快得多。CCG的直观理解是每一次迭代都把最坏的那个分布“定格”下来然后要求主问题在这个固定分布下优化。由于主问题拿到的割是原问题的下界而子问题给出的期望成本是上界上下界不断逼近直到gap小于容差。对一个场景数S15的中等规模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;