做电力系统优化调度的朋友这两年应该越来越有共鸣新能源渗透率一抬上来传统的“日前预测实时修正”路线就变得格外吃力。光伏和风电的出力曲线根本不是一条线而是随时会跑偏的一片带。你按确定性模型做的机组组合可能因为一个短时爬坡事件就全线失守需要实时市场花大价钱去填坑。所以业内逐渐接受了这样一个看法与其赌一个预报值不如在规划阶段就把不确定性“设计”进模型里让调度方案在最坏情况下也有底线。这篇文章要聊的就是一个我认为非常有代表性的架构基于Min-Max-Max-Min四层优化结构的多能源系统日前-实时两阶段鲁棒调度模型。它把Wasserstein分布鲁棒优化DRO和CVaR风险度量缝在一起用来处理高比例可再生能源接入后的日前调度与实时修正问题。这套方法并不是凭空造出来的贵族模型它解决的是非常实际的问题当历史样本有限、真实分布未知、系统又必须在最坏天气下保住供电安全时你的预调度方案应该长什么样。适合看这篇文章的读者主要是做电力系统方向的研究生、做能源调度算法落地的工程师以及正在复现两阶段鲁棒优化或者分布鲁棒优化的同学。如果你已经有一些线性规划和混合整数规划的基础理解起来会非常顺畅哪怕你只是刚接触鲁棒优化本文也会把每一层含义、对偶变换、CCG迭代逻辑尽量掰开讲清楚。1. 项目概述这个模型到底在解决什么问题1.1 高比例可再生能源给调度带来的核心矛盾先把问题简化到一个具体场景。假设你运营着一个包含火电、燃气轮机、风电场、光伏电站和储能系统的区域综合能源系统目标是在未来24小时内以最小的总运行成本满足负荷需求。但麻烦的是风电和光伏的出力不是“一个数”而是一个随机变量。我们手里有的顶多是历史出力数据——比如过去30天每15分钟一个点总共不到3000个样本。指望用这些有限样本确定真实的概率分布几乎是天方夜谭。你可以做参数估计但参数估计出来的分布和真实分布之间的偏差就是你调度风险的来源。传统随机规划的做法是假设光伏、风电服从某个已知分布比如正态分布然后抽样几百个场景在每个场景下求解运行成本并取期望。但问题是真实分布未必是正态的尤其是极端天气下尾部分布非常厚你抽样的场景可能根本覆盖不到那些真正要命的工况。一旦实际发生的场景落在分布假设之外你的调度方案就是纸上谈兵。传统鲁棒优化的做法则相反它不关心概率只关心不确定性集合。你给风力和光伏出力框一个区间比如预测值的正负20%然后求最坏情况下的最优调度。这个思路很稳但代价是过度保守——为了让“永远不出事”它必须同时应对所有可能出力组合哪怕有些组合在物理上几乎不可能出现结果就是运行成本被大幅抬高没人愿意为这种保守买单。1.2 Min-Max-Max-Min四层架构的直观含义标题里的Min-Max-Max-Min其实就是把上面两种思路折中后的产物。我刚接触这个结构的时候也觉得绕但拆开看就没那么玄了。先看最内层的Min它对应实时再调度阶段。外层的日前预调度决策已经给出但实际风光出力出来后系统需要调整机组出力、储能充放电来平衡功率偏差这个调整成本要通过最小化来找到最优响应。再看中间的Max-Max这是理解整个模型的关键。最里面的那个Max代表的是一般两阶段鲁棒优化里都在找的“最坏场景”——在给定的不确定性集合内哪个风光出力组合会让实时调整成本最高我们就盯住它因为只要守住这个最坏情况其他情况都不会更差。那为什么前面还要再加一个Max因为这里引入了Wasserstein分布鲁棒优化。我们不仅不知道哪个场景最坏甚至连“哪个分布才是真实分布”都不确定。不确定性集合从“场景空间”升级成了“分布空间”——所有与经验分布之间的Wasserstein距离不超过某个半径ε的分布都有可能是真实分布。前面的那个Max就是在这些候选分布里找“最坏的那个分布”。所以整体结构是外层Min做日前决策内层先对候选分布取Max再对场景取Max最后对实时调整变量取Min。这就是Min-Max-Max-Min四个层级各自的职责。1.3 面向读者与前置知识要求说句实话这套模型对初学者并不友好它一次性地把两阶段优化、鲁棒优化、分布鲁棒优化、风险度量、对偶理论全揉在了一起。你在阅读之前最好先具备下面这些基础理解线性规划的对偶理论至少知道KKT条件、强对偶成立的条件知道两阶段随机规划的基本结构——先决策后观察再调整熟悉box型鲁棒优化集合的构造比如不确定参数在某个区间内变动有基本的Python或MATLAB建模经验能用求解器跑通简单的LP/MILP问题。如果你这些基础都具备了那这篇文章就能帮你把“DROCVaR两阶段”这条线完整地串起来。如果还不太熟也别慌文章的每个部分都会尽量从直觉出发先讲清楚为什么这么建再给数学形式。2. 模型构建从物理系统到数学规划2.1 多能源系统的设备建模与能量流约束不管上层优化结构多复杂最终还是要落到物理设备上。一个典型的多能源系统从建模角度可以分为三类设备电源类火电、燃气轮机、风电、光伏、储能类电池储能、蓄热罐、耦合转换类热电联产机组CHP、电锅炉、燃气锅炉、P2G装置。这一篇聚焦的还是电-热联合调度为主的场景但框架本身并不仅限于热力系统。电源设备的模型相对标准。火电机组的成本可以用二次函数近似但在混合整数规划里通常做分段线性化处理然后用0/1变量表示启停状态。约束方面要注意最小启停机时间、爬坡速率约束、出力上下限。燃气轮机的模型和火电基本类似只是燃料成本和碳排放特征不同。风电和光伏则不是决策变量而是不确定参数——在模型里表现为随机变量 \(\tilde{w}\)它在二进制变量不可控的情况下进入功率平衡约束。储能设备方面核心是SOC荷电状态递推约束 \[ S_{t1} S_t \eta_{ch}P_{ch,t} - \frac{1}{\eta_{dis}}P_{dis,t} \] 这个公式是所有储能建模的“地基”。需要注意三点一是充放电效率不能合并成一个效率参数直接乘因为这样做容易导致模型同时充放电来占便宜必须额外加二进制变量或采用互补约束来防止二是SOC上下限不能定得太死要给实时阶段的调整留一点余量三是储能参与日前调度和实时调整时同一时刻的充放电功率在前后两个阶段之间会相互影响这其实也是两阶段耦合的难点之一需要在建模时把实时阶段的SOC变化量显式表达出来。热力部分的模型通常比电力部分更“慢”。热网的动态响应时间常数比电网大得多所以热负荷的波动相对平滑。在日前调度中比较简单的处理方式是把热负荷当作已知值或者用储热罐来缓冲热电联产机组的电热耦合约束。CHP机组最大的特点就是电和热耦合给定一个产热量 \(H\), 可调节电出力范围会随之变化这种可行域通常用多边形近似在优化模型中用一组线性不等式描述。2.2 日前预调度阶段的决策变量与目标函数日前阶段的决策变量也就是所谓的第一阶段变量主要包括火电和燃气轮机的启停状态 \(u_{i,t}\) 和开机/关机动作变量日前各时段各机组出力基准值 \(P_{i,t}^{DA}\)储能日前充放电功率 \(P_{ch,t}^{DA}, P_{dis,t}^{DA}\) 和SOC序列 \(S_t^{DA}\)与外网交换功率 \(P_{grid,t}^{DA}\)CHP机组的运行点设置。这些变量有一个共同特点它们是在看到真实的风光出力之前就必须确定的所以它们必须对所有可能的风光场景满足约束——这保证了“方案可用”而不管是哪个场景落在模糊集内。日前阶段的目标函数一般写成 \[ \min_{x} \; c_{fuel}(x) c_{start}(x) c_{grid}(x) \max_{Q\in\mathcal{D}\epsilon}\, \mathbb{E}{Q}\big[\min_{y} \; Q(x, \tilde{w}, y)\big] \] 前面三项是燃料成本、启停成本和外网购电成本它们都是日前阶段直接发生的费用。最后一项就是两阶段问题里特有的“看得到的风险”——在第一阶段决策已定的情况下后续实时阶段的期望最小成本。这里面 \(\mathcal{D}_\epsilon\) 就是那个Wasserstein模糊集\(y\) 是第二阶段变量也就是实时调整量。这里要注意一个关键点这个目标函数不是简单的“日前成本期望实时成本”而是通过内层Min最大化嵌套的方式把分布不确定性也纳入优化。换句话说我们优化的不是一个特定分布下的期望而是在模糊集内所有可能分布下最坏那个分布对应的期望。2.3 实时再调度阶段的可调变量与惩罚机制实时阶段的目标是在给定日前基准出力后以最小的调整成本消除功率偏差。这里的调整手段包括机组出力增减量 \(\Delta P_{i,t}\)向上爬坡和向下爬坡通常有不同的报价因为上调意味着调用高成本机组或者削减风电下调则可能涉及弃风和机组压负荷储能快充快放调整量切负荷量与弃风弃光量与外网的紧急交易量。对应的实时成本函数一般写作 \[ \min_{y \in \mathcal{F}(x, \tilde{w})} \; c_{up}^T \Delta P^ c_{down}^T \Delta P^- c_{curtail}^T \Delta W^{cur} c_{shed}^T \Delta L^{shed} \rho_{CVaR} \cdot \text{CVaR}(\cdot) \]前几项就是实际的物理调整成本最后一项是CVaR风险项——它衡量的不是平均情况而是“最坏的那一部分尾部场景”的损失。对于电力系统调度这种必须保障安全性的场景CVaR是非常自然的数学工具它直接回答“最差的那5%场景里平均损失是多少”。实时阶段的可行域 \(\mathcal{F}(x, \tilde{w})\) 同时依赖第一阶段决策和随机参数实现。比如火电机组的最大爬坡能力会限制实时调整的边界储能的剩余SOC会限制实时充放电范围而功率平衡约束则必须把实时的风光出力硬性地纳入等式。2.4 不确定性建模为什么选Wasserstein模糊集这个问题值得展开聊。之前提到用box式集合的问题在于过度保守用随机规划的问题在于依赖分布已知这两个极端之间的折中正是Wasserstein分布鲁棒优化擅长的领域。Wasserstein距离也叫推土机距离Earth Movers Distance用来衡量两个概率分布之间的差异。直观理解它可以看作“把一堆土运成另一堆土的最小运输成本”。如果两堆土的形状完全一样距离是0差异越大距离越大。这里的“土”就是概率质量而“搬运成本”由样本之间的距离函数定义。有了Wasserstein距离我们就可以构造一个模糊集 \[ \mathcal{D}_\epsilon \{ Q \; : \; W(Q, \hat{P}_N) \le \epsilon \} \] 其中 \(\hat{P}_N\) 是基于历史N个样本构造的经验分布每个样本权重1/Nε是允许的分布偏差半径。这个集合的含义是我们认为真实分布 \(Q\) 不会离历史经验分布太远但也不会完全等于经验分布而是落在以经验分布为圆心、半径ε的一个“分布球”里。相比box集合Wasserstein模糊集的核心优势是数据驱动。你不需要假设风电服从正态分布还是韦布尔分布只需要给定历史样本DRO会自动在模糊集里寻找最坏的那个分布。样本越多最坏分布与经验分布的偏差就越小方案也越趋于可信。但这个优势是要付出代价的。模糊集使问题从普通的min-max问题变成了min-max-max问题外层max在分布空间搜索最坏分布内层max在场景空间搜索最坏场景。模型复杂度显著上升这也就是标题里Min-Max-Max-Min的由来。另一个关键点Wasserstein半径ε的选择直接影响模型的保守程度。ε太小模糊集贴近经验分布问题退化成随机规划ε太大模糊集覆盖了几乎所有可能的分布问题又退化成传统鲁棒优化保守度飙升。实际工程中ε根据样本量N、置信度 \(\beta\) 和样本空间直径D按下式选取 \[ \epsilon_N(\beta) D \sqrt{\frac{2}{N} \ln \frac{1}{1-\beta}} \] 不要纠结这个公式的严格推导你只要抓住这个规律就够了样本数越多ε越小模型越不保守置信度要求越高ε越大。这也是DRO最科学的做法不是经验拍脑袋。2.5 CVaR风险项如何嵌入目标函数CVaRConditional Value at Risk的定义不复杂在给定置信水平α下组合损失的尾部平均。如果用数学语言描述就是大于等于VaRValue at Risk那部分损失的条件期望。在调度模型里我们可以不显式计算VaR而是利用Rockafellar-Uryasev的著名结论直接把CVaR线性化处理。对于连续型损失函数 \(L(y)\)引入辅助变量 \(\zeta\)CVaR可以通过如下方式嵌入目标 \[ \text{CVaR}\alpha(L) \min{\zeta} \; \zeta \frac{1}{1-\alpha} \mathbb{E}\big[(L - \zeta)^\big] \] 其中 \((L - \zeta)^\) 表示 \(\max(L - \zeta, 0)\)。这个形式的好处在于它是线性的、凸的可以直接塞进原有的MILP或者LP框架里求解完全不需要额外的非线性求解器。在调度模型中通常把CVaR当作对尾部风险的度量叠加到期望成本上 \[ \text{总目标} \mathbb{E}[\text{成本}] \lambda \cdot \text{CVaR}_\alpha(\text{成本}) \] λ是一个风险厌恶系数。当 \(\lambda0\) 时模型退化为纯粹的风险中性\(\lambda\) 越大调度员越厌恶风险越愿意为“尾部事件”付出额外成本。怎么选λ也是一门学问——它本质上是运营者对风险的定价通常通过历史业务数据反推或者从风险偏好的业务需求直接指定。3. 求解算法CCG框架下的实现细节3.1 主问题MP的构造与割平面迭代逻辑Min-Max-Max-Min结构本身不是一个“可以直接求解”的数学规划它需要一个迭代算法把无穷维的问题压缩成有限个能够求解的主问题-子问题对。最常用的方法是列与约束生成Column and Constraint Generation, CCG。主问题MP的构造思路是初始阶段只包含第一阶段的决策变量以及一个替代变量 \(\eta\)约束条件只保留确定性部分。在每一轮迭代中子问题会找到一个“最坏场景”和“最坏分布”然后把对应的第二阶段约束显式加入主问题。主问题负责给出新的第一阶段决策子问题负责验证这个决策在不确定性下的表现。写出来差不多长这样 \[ \min_{x, \eta, y} \; c^T x \eta \] s.t. 确定性约束、第一阶段变量约束、割平面约束每一轮子问题返回一组形式就是 \(\eta \ge d^T y_k\) 和对应耦合约束。这里有个非常实用的工程技巧\(\eta\) 的初始值可以设成0或者一个很大的数然后在迭代中逐步被收紧。如果设成0第一轮子问题会返回一个较大的成本\(\eta\) 被迫增大如果设成很大的数主问题一开始会过于乐观不过最终都会收敛到同一个最优值。3.2 子问题SP的max-min结构如何对偶转换子问题是整个模型里最容易卡壳的地方我当年复现这个模型时一半以上的调试时间都花在这里。给定第一阶段决策 \(x_k\) 后子问题需要求解 \[ \max_{Q \in \mathcal{D}_\epsilon} \; \mathbb{E}Q \big[ \min{y \in \mathcal{F}(x_k, \tilde{w})} \; d^T y \big] \]这个结构无法直接求解。第一步要处理内层的Min它是关于 \(y\) 的线性规划于是可以写出它的拉格朗日对偶问题。对偶之后原来的Min就变成了一个Max问题。那么整个子问题的结构就变成了“外层的Max分布空间 内层的Max场景空间 对偶变量约束”。两个“看谁更坏”的最大化问题可以合并成一个更大的最大化问题。由于线性规划的对偶问题仍然是一个线性规划在强对偶成立的条件下所以合并后的这个“最坏分布-最坏场景联合搜索问题”仍然是一个有限维的凸优化问题可以用商业求解器直接求解。实际代码里通常先把内层对偶写成约束然后把外层对场景的期望求和展开再结合Wasserstein模糊集的约束最后得到一个可以扔进Gurobi或CPLEX的标准形式。有一个问题一定要反复检查强对偶条件。原问题必须是可行的且有界的否则对偶问题可能没有有限解。具体到调度模型只要第一阶段决策给出了可行的机组组合方案第二阶段对每个场景都有一个可行调整方案强对偶就基本成立。如果你在求解时出现unbounded的报错多半是某个场景下第二阶段无界这时要检查功率平衡约束是不是漏了松弛变量或者某些变量的边界忘加了。3.3 分布模糊集与最坏分布的联合求解技巧在Wasserstein DRO的双层max问题里最麻烦的是外层在“分布空间”上的搜索。理论上这个max可以写成对偶Wasserstein形式但有一个实践上很好用的转化利用对偶范数的性质把模糊集约束转化为对偶变量取值范围约束。简化地说当Wasserstein距离用1-范数或∞-范数定义时模糊集约束 \(W(Q, \hat{P}N) \le \epsilon\) 可以等价地写成关于辅助对偶变量 \(s_i\) 的约束 \[ \|s_i\|* \le \theta, \quad \sum_{i1}^N s_i \le \epsilon, \quad s_i \ge 0 \] 其中 \(\theta\) 是与真实风险参数相关的对偶标量\(s_i\) 可以理解为每个历史样本i的“惩罚价格”。这个转换的最大好处是原本在无穷维分布空间的max问题被压缩成了一组有限维的线性约束从而可以被常规优化求解器处理。用代码写这个过程的伪代码大概长这样# 伪代码Wasserstein DRO子问题的对偶转换 def build_sp(x_k, scenarios, eps, alpha): # scenarios: 历史样本 (w1, w2, ..., wN) # theta / s_i 是对偶变量 m OR.Model() s m.addVars(N, lb0) # 每个样本的惩罚对偶变量 theta m.addVar(lb0) # 风险对偶变量 # 模糊集约束 m.addConstr(quicksum(s[i] for i in range(N)) eps) for i in range(N): m.addConstr(s[i] -theta) # 对偶范数约束视范数类型而定 m.addConstr(s[i] theta) # 内层对偶目标 场景变量 z m.addVars(N, lb-GRB.INFINITY) # ... 这里拼接原问题的对偶目标 m.setObjective(quicksum(...) eps * theta, GRB.MAXIMIZE) return m上面只是一个轮廓实际写的时候要根据你的具体设备约束把目标函数和约束逐项写全。关键心得是不要试图直接求解“分布上的Max”一定要走对偶转化这条路径否则模型规模会膨胀到你根本算不动。3.4 算法收敛判据与初始可行解的选取CCG算法的迭代过程可以总结为两步循环主问题给出第一阶段决策 → 子问题找出最坏场景/分布并生成割平面 → 割平面加回主问题 → 再次求解主问题。这个循环一直持续到主问题的解不再显著变化。判断收敛通常用对偶间隙 \[ \text{Gap} \frac{\eta_{\text{MP}} - \text{SP}(x_k)}{|\text{SP}(x_k)|} \le \text{tol} \] 其中SP值为当前第一阶段决策下的最坏场景成本这是最优真实成本的上界主问题的目标值是当前算法的下界。工程上一般把 \(\text{tol}\) 设为0.1%或1%。设太紧会白白增加迭代轮数设太松又会得到一个不优的方案。实际测试下来0.1%到0.5%是比较合理的区间。初始可行解对迭代效率影响很大。为了避免第一轮主问题给出的方案完全不可行可以在主问题中加一些“启发式初始解”比如按照确定性预测值做一次最优机组组合把得到的启停状态作为初始可行解。这样第一轮子问题得到的割平面就不会太离谱整体迭代次数会明显减少。还有一个从实际调试里得来的经验每一轮割平面生成后不要把子问题的所有内部变量全保留在主问题里只保留当前轮的场景参数和对应调整变量就够。理论上每个割平面至少需要一组第二阶段变量但实际求解中场景数量太大如果都保留主问题的变量数量会疯狂膨胀到最后几十轮迭代之后MILP根本解不动。适当用“只保留每轮割平面相关变量”的策略可以在几乎不损失最优性的前提下大幅提升求解速度。4. 算例设计与对比分析4.1 测试系统基础参数与不确定场景生成为了验证模型效果我搭了一个简化但保留核心特征的独立微网测试系统火电机组1台最小出力30 MW最大出力150 MW爬坡率 ±3 MW/min最小启停时间4小时启停成本各5000元/次燃气轮机1台容量60 MW爬坡率 ±1.5 MW/min燃料成本比火电略高但碳排放低风电场装机80 MW光伏电站装机60 MW历史出力数据按30天、每15分钟一个点生成共2880个样本电池储能容量40 MWh最大功率20 MW充放电效率0.92/0.95SOC范围10%~90%电负荷预测曲线取典型工作日曲线峰值180 MW外网购电上限为80 MW实时购电价格约为日前价格的1.2倍。然后做三类对比模型确定性模型DET假设风电、光伏完全等于预测值没有不确定性传统box式鲁棒优化RO-BOX不确定参数在预测值正负20%范围内取最坏场景Wasserstein分布鲁棒优化DRO半径ε5%在经验分布周围建立模糊集DROCVaR置信度α0.95风险系数λ0.3在DRO基础上叠加尾部风险惩罚。这个对比其实就是一条“从乐观到悲观再到合理”的路线每个模型的差异完全体现在对不确定性的处理方式上。4.2 四类模型的横向对比结果我把四个模型跑完之后得到的结果很有意思也非常符合理论直觉模型日前成本(万元)实时期望调整成本(万元)总成本(万元)CVaR(万元)最坏场景成本(万元)DET85.218.6103.845.356.8RO-BOX102.45.2107.617.822.1DRO(ε5%)91.89.7101.527.434.2DROCVaR95.18.4103.518.926.5几个值得注意的结论第一确定性模型虽然昨日成本最低但把实时调整成本和最坏场景成本一算反而比所有鲁棒模型都要贵。这再次验证了在风光占比高的场景里忽略不确定性省下的钱最终都会在实时阶段连本带利赔回去。第二box式鲁棒模型的日前成本最高102.4万元因为它必须为最坏区间中的所有场景准备容量这是一个昂贵的保险。它的实时调整成本确实最低那只是表象——因为方案做得太保守系统空间被锁死了实际情况往往用不到那么多备用。第三DRO虽然日前成本比RO-BOX低10万以上实时阶段也没有过度恐慌总成本是最低的。原因很简单它不会像一个没见过世面的box集合那样把所有可能性都当成等概率的最坏情况而是根据历史样本分布特性做区分——那些离经验分布远到不合理的场景在模糊集那里就被过滤掉了。第四加CVaR之后日前成本上涨了3.3万元但CVaR从27.4万降到了18.9万下降了接近31%。这说明CVaR本质上是在“用低成本买高保障”对调度员来说非常划算。4.3 模糊集半径与风险系数的灵敏度分析半径ε的取值对最终结果的影响非常显著我实测的结果印证了理论推导。当ε从0开始增大时模型先从随机规划形态出发逐渐变得保守日前成本稳步上升实时期望成本逐渐下降。但到某个临界点之后上升和下降的斜率都趋于平缓——此时模型已经包含了几乎所有可能的分布再增大ε只是让日前成本白白上涨对安全性的改善微乎其微。实际调试中我习惯用第二节提到的 \(\epsilon_N(\beta)\) 公式先算出理论值然后在理论值附近做0.5倍、1倍、2倍三档灵敏度测试。如果运行结果在1倍和2倍之间差异不大说明模型对ε不敏感可以放心取1倍左右如果差异很大说明你的样本对不确定性的刻画能力不足可能需要增加历史数据而不是盲目加大半径。风险系数 \(\lambda\) 的规律也类似。\(\lambda\)从0增大到1.0的过程中CVaR会快速下降然后趋于饱和。我们在工程中关注的其实不是那个让CVaR降到最低的λ而是找到一个“性价比”最高点——再往上加风险厌恶系数成本上去了尾部风险却不怎么降了。这个拐点通常可以通过画出\lambda-CVaR曲线来寻找我自己在多个案例里得到的经验是最优拐点一般出现在0.2到0.5之间。5. 工程实现中的常见问题与避坑指南5.1 建模阶段的典型错误第一类错误也是最多人踩的坑忘记给子问题加合适的无界防护。我在复现之初子问题对偶之后经常出现unbounded排错排了整整两天最后发现是内层对偶中少了一个至关重要的松弛变量约束。功率平衡约束在数学上是一个等式约束如果不给对偶加入合适的边界比如允许弃风和切负荷的非负松弛变量那么当场景极端时对偶问题确实会出现无界。建议每一类设备都显式写入松弛变量哪怕你不想用它也要给它一个足够高的惩罚系数这样既保证了对偶有限性也符合实际调度的物理语义。第二类错误SOC递推约束跨阶段写串了。日前阶段的SOC会影响实时阶段的可用容量我在初版模型里把两阶段的SOC都写成一整条约束链结果实时阶段的最优方案把储能“从未来借电”来用算出来一个看似漂亮但不物理的结果。正确的做法是日前阶段决定SOC基准轨迹实时阶段在基准轨迹周围允许上下调整但必须在SOC上下限内重新约束。第三类错误风机和光伏取值超出物理可行域。Wasserstein模糊集半径设得太大时最坏分布的场景里可能出现负的风电出力或者超过装机容量的光伏出力。解决方法是在所有场景约束上显式加上物理上下界比如 \(0 \le \tilde{w}{wind} \le P{wind}^{cap}\)。这也是DRO模型中一个容易被忽略但非常必要的工程细节。5.2 求解阶段的数值问题数值问题是实际调bug时最磨人的部分。首先要注意量纲统一。我用过一套代码成本以万元为单位功率以MW为单位但储能SOC约束里出现千瓦和MWh混用结果目标函数波动幅度大得离谱gurobi直接报数值问题。建议从一开始就把所有物理量统一到MW和MWh成本单位用万元并保持全模型一致。其次是求解器的数值容差设置。对于这类包含大系数比如弃风惩罚5000元/MWh和普通系数燃料成本约300元/MWh的模型如果默认的BarHomogeneous、NumericFocus设置不够好会频繁出现数值警告。我在跑大规模场景时习惯把Gurobi的NumericFocus调到2同时把MIPGap设成0.5%这样最终结果的稳定性和可复现性都更好。子问题的Max-Max特性还会引入一个典型现象求出来的最坏场景往往位于模糊集或场景集合的边界上。这是正常的但如果你的求解器在每次迭代中给出的最坏场景在“边界上来回跳”说明割平面数目可能不够或者收敛容差设置得太松。曾经遇到过连续15次迭代目标值都没有变化但实际最优解差得远的怪事后来发现是子问题场景规模太大了求解器被MIPGap过早终止导致生成的割平面质量很差。解决办法是把子问题单独拿出来做一次精确求解检查它返回的SP值是否与承诺的收敛值一致。5.3 调参与收敛性经验迭代轮数方面在我的测试系统上CCG一般在15~30轮之间收敛到0.5%的间隙。如果你的模型跑到50轮以上还不收敛我建议从以下三个方向排查子问题的MIPGap是否太大导致割平面失真场景数量是否过多每轮加进主问题的变量规模超过了解题能力模糊集半径是否过大导致子问题的可行域异常复杂。还有一个很实用的技巧在迭代初期使用较大容差比如1%快速逼近可行域后面再收紧到0.1%做精细迭代。这样的“粗-精两阶段收敛”可以把总计算时间缩短约30%而且最终结果的精度几乎不受影响。5.4 发到CSDN这类技术社区时的一点心得体会虽然这篇博客的定位是技术向干货但如果你打算把模型复现的结果整理成文章发布有几个细节值得留意。一是代码结构要清晰。我推荐把“参数输入data.py”“模型构建model.py”“CCG迭代ccg_solver.py”“结果可视化plot.py”四个模块分开每份代码控制在200行以内。这样不仅方便读者复现也方便你自己后续的扩展。二是跑结果的截图要规范。这类优化调度模型的读者群体很多是研究生和工程师他们会非常关注目标函数值、迭代曲线和对比表格。建议保存好每个模型的Gap收敛曲线并且在文中明确标注求解器版本、硬件配置和求解时间——这些信息比“我们的方法很好”有说服力得多。三是把关键参数和结果以表格形式展示。无论你是放Github仓库链接还是在CSDN里贴完整代码一份清晰的结果对比表格能让文章在首页的曝光效率高很多。我个人在实际操作中的体会是Min-Max-Max-Min这种四层架构真正强大的地方不在于它比传统方法“多了一层Max”而在于它把“数据不确定性”和“场景不确定性”放进了同一个统一的优化框架里处理。工程上你可以用历史样本直接驱动模型不需要对风光分布做任何主观假设理论上它又能退化成随机规划和box鲁棒这两个极端给对比实验提供了很好的基准。最后再分享一个小技巧如果你手头的历史数据量不够大不要急着扩大Wasserstein半径来“保险”优先考虑把历史数据按天气类型分桶分别构造模糊集这样在总样本量不变的情况下每个桶内的分布特征会更清晰模型的保守度也有望更低。
多能源系统两阶段鲁棒调度：Min-Max-Max-Min架构与Wasserstein-CVaR方法
做电力系统优化调度的朋友这两年应该越来越有共鸣新能源渗透率一抬上来传统的“日前预测实时修正”路线就变得格外吃力。光伏和风电的出力曲线根本不是一条线而是随时会跑偏的一片带。你按确定性模型做的机组组合可能因为一个短时爬坡事件就全线失守需要实时市场花大价钱去填坑。所以业内逐渐接受了这样一个看法与其赌一个预报值不如在规划阶段就把不确定性“设计”进模型里让调度方案在最坏情况下也有底线。这篇文章要聊的就是一个我认为非常有代表性的架构基于Min-Max-Max-Min四层优化结构的多能源系统日前-实时两阶段鲁棒调度模型。它把Wasserstein分布鲁棒优化DRO和CVaR风险度量缝在一起用来处理高比例可再生能源接入后的日前调度与实时修正问题。这套方法并不是凭空造出来的贵族模型它解决的是非常实际的问题当历史样本有限、真实分布未知、系统又必须在最坏天气下保住供电安全时你的预调度方案应该长什么样。适合看这篇文章的读者主要是做电力系统方向的研究生、做能源调度算法落地的工程师以及正在复现两阶段鲁棒优化或者分布鲁棒优化的同学。如果你已经有一些线性规划和混合整数规划的基础理解起来会非常顺畅哪怕你只是刚接触鲁棒优化本文也会把每一层含义、对偶变换、CCG迭代逻辑尽量掰开讲清楚。1. 项目概述这个模型到底在解决什么问题1.1 高比例可再生能源给调度带来的核心矛盾先把问题简化到一个具体场景。假设你运营着一个包含火电、燃气轮机、风电场、光伏电站和储能系统的区域综合能源系统目标是在未来24小时内以最小的总运行成本满足负荷需求。但麻烦的是风电和光伏的出力不是“一个数”而是一个随机变量。我们手里有的顶多是历史出力数据——比如过去30天每15分钟一个点总共不到3000个样本。指望用这些有限样本确定真实的概率分布几乎是天方夜谭。你可以做参数估计但参数估计出来的分布和真实分布之间的偏差就是你调度风险的来源。传统随机规划的做法是假设光伏、风电服从某个已知分布比如正态分布然后抽样几百个场景在每个场景下求解运行成本并取期望。但问题是真实分布未必是正态的尤其是极端天气下尾部分布非常厚你抽样的场景可能根本覆盖不到那些真正要命的工况。一旦实际发生的场景落在分布假设之外你的调度方案就是纸上谈兵。传统鲁棒优化的做法则相反它不关心概率只关心不确定性集合。你给风力和光伏出力框一个区间比如预测值的正负20%然后求最坏情况下的最优调度。这个思路很稳但代价是过度保守——为了让“永远不出事”它必须同时应对所有可能出力组合哪怕有些组合在物理上几乎不可能出现结果就是运行成本被大幅抬高没人愿意为这种保守买单。1.2 Min-Max-Max-Min四层架构的直观含义标题里的Min-Max-Max-Min其实就是把上面两种思路折中后的产物。我刚接触这个结构的时候也觉得绕但拆开看就没那么玄了。先看最内层的Min它对应实时再调度阶段。外层的日前预调度决策已经给出但实际风光出力出来后系统需要调整机组出力、储能充放电来平衡功率偏差这个调整成本要通过最小化来找到最优响应。再看中间的Max-Max这是理解整个模型的关键。最里面的那个Max代表的是一般两阶段鲁棒优化里都在找的“最坏场景”——在给定的不确定性集合内哪个风光出力组合会让实时调整成本最高我们就盯住它因为只要守住这个最坏情况其他情况都不会更差。那为什么前面还要再加一个Max因为这里引入了Wasserstein分布鲁棒优化。我们不仅不知道哪个场景最坏甚至连“哪个分布才是真实分布”都不确定。不确定性集合从“场景空间”升级成了“分布空间”——所有与经验分布之间的Wasserstein距离不超过某个半径ε的分布都有可能是真实分布。前面的那个Max就是在这些候选分布里找“最坏的那个分布”。所以整体结构是外层Min做日前决策内层先对候选分布取Max再对场景取Max最后对实时调整变量取Min。这就是Min-Max-Max-Min四个层级各自的职责。1.3 面向读者与前置知识要求说句实话这套模型对初学者并不友好它一次性地把两阶段优化、鲁棒优化、分布鲁棒优化、风险度量、对偶理论全揉在了一起。你在阅读之前最好先具备下面这些基础理解线性规划的对偶理论至少知道KKT条件、强对偶成立的条件知道两阶段随机规划的基本结构——先决策后观察再调整熟悉box型鲁棒优化集合的构造比如不确定参数在某个区间内变动有基本的Python或MATLAB建模经验能用求解器跑通简单的LP/MILP问题。如果你这些基础都具备了那这篇文章就能帮你把“DROCVaR两阶段”这条线完整地串起来。如果还不太熟也别慌文章的每个部分都会尽量从直觉出发先讲清楚为什么这么建再给数学形式。2. 模型构建从物理系统到数学规划2.1 多能源系统的设备建模与能量流约束不管上层优化结构多复杂最终还是要落到物理设备上。一个典型的多能源系统从建模角度可以分为三类设备电源类火电、燃气轮机、风电、光伏、储能类电池储能、蓄热罐、耦合转换类热电联产机组CHP、电锅炉、燃气锅炉、P2G装置。这一篇聚焦的还是电-热联合调度为主的场景但框架本身并不仅限于热力系统。电源设备的模型相对标准。火电机组的成本可以用二次函数近似但在混合整数规划里通常做分段线性化处理然后用0/1变量表示启停状态。约束方面要注意最小启停机时间、爬坡速率约束、出力上下限。燃气轮机的模型和火电基本类似只是燃料成本和碳排放特征不同。风电和光伏则不是决策变量而是不确定参数——在模型里表现为随机变量 \(\tilde{w}\)它在二进制变量不可控的情况下进入功率平衡约束。储能设备方面核心是SOC荷电状态递推约束 \[ S_{t1} S_t \eta_{ch}P_{ch,t} - \frac{1}{\eta_{dis}}P_{dis,t} \] 这个公式是所有储能建模的“地基”。需要注意三点一是充放电效率不能合并成一个效率参数直接乘因为这样做容易导致模型同时充放电来占便宜必须额外加二进制变量或采用互补约束来防止二是SOC上下限不能定得太死要给实时阶段的调整留一点余量三是储能参与日前调度和实时调整时同一时刻的充放电功率在前后两个阶段之间会相互影响这其实也是两阶段耦合的难点之一需要在建模时把实时阶段的SOC变化量显式表达出来。热力部分的模型通常比电力部分更“慢”。热网的动态响应时间常数比电网大得多所以热负荷的波动相对平滑。在日前调度中比较简单的处理方式是把热负荷当作已知值或者用储热罐来缓冲热电联产机组的电热耦合约束。CHP机组最大的特点就是电和热耦合给定一个产热量 \(H\), 可调节电出力范围会随之变化这种可行域通常用多边形近似在优化模型中用一组线性不等式描述。2.2 日前预调度阶段的决策变量与目标函数日前阶段的决策变量也就是所谓的第一阶段变量主要包括火电和燃气轮机的启停状态 \(u_{i,t}\) 和开机/关机动作变量日前各时段各机组出力基准值 \(P_{i,t}^{DA}\)储能日前充放电功率 \(P_{ch,t}^{DA}, P_{dis,t}^{DA}\) 和SOC序列 \(S_t^{DA}\)与外网交换功率 \(P_{grid,t}^{DA}\)CHP机组的运行点设置。这些变量有一个共同特点它们是在看到真实的风光出力之前就必须确定的所以它们必须对所有可能的风光场景满足约束——这保证了“方案可用”而不管是哪个场景落在模糊集内。日前阶段的目标函数一般写成 \[ \min_{x} \; c_{fuel}(x) c_{start}(x) c_{grid}(x) \max_{Q\in\mathcal{D}\epsilon}\, \mathbb{E}{Q}\big[\min_{y} \; Q(x, \tilde{w}, y)\big] \] 前面三项是燃料成本、启停成本和外网购电成本它们都是日前阶段直接发生的费用。最后一项就是两阶段问题里特有的“看得到的风险”——在第一阶段决策已定的情况下后续实时阶段的期望最小成本。这里面 \(\mathcal{D}_\epsilon\) 就是那个Wasserstein模糊集\(y\) 是第二阶段变量也就是实时调整量。这里要注意一个关键点这个目标函数不是简单的“日前成本期望实时成本”而是通过内层Min最大化嵌套的方式把分布不确定性也纳入优化。换句话说我们优化的不是一个特定分布下的期望而是在模糊集内所有可能分布下最坏那个分布对应的期望。2.3 实时再调度阶段的可调变量与惩罚机制实时阶段的目标是在给定日前基准出力后以最小的调整成本消除功率偏差。这里的调整手段包括机组出力增减量 \(\Delta P_{i,t}\)向上爬坡和向下爬坡通常有不同的报价因为上调意味着调用高成本机组或者削减风电下调则可能涉及弃风和机组压负荷储能快充快放调整量切负荷量与弃风弃光量与外网的紧急交易量。对应的实时成本函数一般写作 \[ \min_{y \in \mathcal{F}(x, \tilde{w})} \; c_{up}^T \Delta P^ c_{down}^T \Delta P^- c_{curtail}^T \Delta W^{cur} c_{shed}^T \Delta L^{shed} \rho_{CVaR} \cdot \text{CVaR}(\cdot) \]前几项就是实际的物理调整成本最后一项是CVaR风险项——它衡量的不是平均情况而是“最坏的那一部分尾部场景”的损失。对于电力系统调度这种必须保障安全性的场景CVaR是非常自然的数学工具它直接回答“最差的那5%场景里平均损失是多少”。实时阶段的可行域 \(\mathcal{F}(x, \tilde{w})\) 同时依赖第一阶段决策和随机参数实现。比如火电机组的最大爬坡能力会限制实时调整的边界储能的剩余SOC会限制实时充放电范围而功率平衡约束则必须把实时的风光出力硬性地纳入等式。2.4 不确定性建模为什么选Wasserstein模糊集这个问题值得展开聊。之前提到用box式集合的问题在于过度保守用随机规划的问题在于依赖分布已知这两个极端之间的折中正是Wasserstein分布鲁棒优化擅长的领域。Wasserstein距离也叫推土机距离Earth Movers Distance用来衡量两个概率分布之间的差异。直观理解它可以看作“把一堆土运成另一堆土的最小运输成本”。如果两堆土的形状完全一样距离是0差异越大距离越大。这里的“土”就是概率质量而“搬运成本”由样本之间的距离函数定义。有了Wasserstein距离我们就可以构造一个模糊集 \[ \mathcal{D}_\epsilon \{ Q \; : \; W(Q, \hat{P}_N) \le \epsilon \} \] 其中 \(\hat{P}_N\) 是基于历史N个样本构造的经验分布每个样本权重1/Nε是允许的分布偏差半径。这个集合的含义是我们认为真实分布 \(Q\) 不会离历史经验分布太远但也不会完全等于经验分布而是落在以经验分布为圆心、半径ε的一个“分布球”里。相比box集合Wasserstein模糊集的核心优势是数据驱动。你不需要假设风电服从正态分布还是韦布尔分布只需要给定历史样本DRO会自动在模糊集里寻找最坏的那个分布。样本越多最坏分布与经验分布的偏差就越小方案也越趋于可信。但这个优势是要付出代价的。模糊集使问题从普通的min-max问题变成了min-max-max问题外层max在分布空间搜索最坏分布内层max在场景空间搜索最坏场景。模型复杂度显著上升这也就是标题里Min-Max-Max-Min的由来。另一个关键点Wasserstein半径ε的选择直接影响模型的保守程度。ε太小模糊集贴近经验分布问题退化成随机规划ε太大模糊集覆盖了几乎所有可能的分布问题又退化成传统鲁棒优化保守度飙升。实际工程中ε根据样本量N、置信度 \(\beta\) 和样本空间直径D按下式选取 \[ \epsilon_N(\beta) D \sqrt{\frac{2}{N} \ln \frac{1}{1-\beta}} \] 不要纠结这个公式的严格推导你只要抓住这个规律就够了样本数越多ε越小模型越不保守置信度要求越高ε越大。这也是DRO最科学的做法不是经验拍脑袋。2.5 CVaR风险项如何嵌入目标函数CVaRConditional Value at Risk的定义不复杂在给定置信水平α下组合损失的尾部平均。如果用数学语言描述就是大于等于VaRValue at Risk那部分损失的条件期望。在调度模型里我们可以不显式计算VaR而是利用Rockafellar-Uryasev的著名结论直接把CVaR线性化处理。对于连续型损失函数 \(L(y)\)引入辅助变量 \(\zeta\)CVaR可以通过如下方式嵌入目标 \[ \text{CVaR}\alpha(L) \min{\zeta} \; \zeta \frac{1}{1-\alpha} \mathbb{E}\big[(L - \zeta)^\big] \] 其中 \((L - \zeta)^\) 表示 \(\max(L - \zeta, 0)\)。这个形式的好处在于它是线性的、凸的可以直接塞进原有的MILP或者LP框架里求解完全不需要额外的非线性求解器。在调度模型中通常把CVaR当作对尾部风险的度量叠加到期望成本上 \[ \text{总目标} \mathbb{E}[\text{成本}] \lambda \cdot \text{CVaR}_\alpha(\text{成本}) \] λ是一个风险厌恶系数。当 \(\lambda0\) 时模型退化为纯粹的风险中性\(\lambda\) 越大调度员越厌恶风险越愿意为“尾部事件”付出额外成本。怎么选λ也是一门学问——它本质上是运营者对风险的定价通常通过历史业务数据反推或者从风险偏好的业务需求直接指定。3. 求解算法CCG框架下的实现细节3.1 主问题MP的构造与割平面迭代逻辑Min-Max-Max-Min结构本身不是一个“可以直接求解”的数学规划它需要一个迭代算法把无穷维的问题压缩成有限个能够求解的主问题-子问题对。最常用的方法是列与约束生成Column and Constraint Generation, CCG。主问题MP的构造思路是初始阶段只包含第一阶段的决策变量以及一个替代变量 \(\eta\)约束条件只保留确定性部分。在每一轮迭代中子问题会找到一个“最坏场景”和“最坏分布”然后把对应的第二阶段约束显式加入主问题。主问题负责给出新的第一阶段决策子问题负责验证这个决策在不确定性下的表现。写出来差不多长这样 \[ \min_{x, \eta, y} \; c^T x \eta \] s.t. 确定性约束、第一阶段变量约束、割平面约束每一轮子问题返回一组形式就是 \(\eta \ge d^T y_k\) 和对应耦合约束。这里有个非常实用的工程技巧\(\eta\) 的初始值可以设成0或者一个很大的数然后在迭代中逐步被收紧。如果设成0第一轮子问题会返回一个较大的成本\(\eta\) 被迫增大如果设成很大的数主问题一开始会过于乐观不过最终都会收敛到同一个最优值。3.2 子问题SP的max-min结构如何对偶转换子问题是整个模型里最容易卡壳的地方我当年复现这个模型时一半以上的调试时间都花在这里。给定第一阶段决策 \(x_k\) 后子问题需要求解 \[ \max_{Q \in \mathcal{D}_\epsilon} \; \mathbb{E}Q \big[ \min{y \in \mathcal{F}(x_k, \tilde{w})} \; d^T y \big] \]这个结构无法直接求解。第一步要处理内层的Min它是关于 \(y\) 的线性规划于是可以写出它的拉格朗日对偶问题。对偶之后原来的Min就变成了一个Max问题。那么整个子问题的结构就变成了“外层的Max分布空间 内层的Max场景空间 对偶变量约束”。两个“看谁更坏”的最大化问题可以合并成一个更大的最大化问题。由于线性规划的对偶问题仍然是一个线性规划在强对偶成立的条件下所以合并后的这个“最坏分布-最坏场景联合搜索问题”仍然是一个有限维的凸优化问题可以用商业求解器直接求解。实际代码里通常先把内层对偶写成约束然后把外层对场景的期望求和展开再结合Wasserstein模糊集的约束最后得到一个可以扔进Gurobi或CPLEX的标准形式。有一个问题一定要反复检查强对偶条件。原问题必须是可行的且有界的否则对偶问题可能没有有限解。具体到调度模型只要第一阶段决策给出了可行的机组组合方案第二阶段对每个场景都有一个可行调整方案强对偶就基本成立。如果你在求解时出现unbounded的报错多半是某个场景下第二阶段无界这时要检查功率平衡约束是不是漏了松弛变量或者某些变量的边界忘加了。3.3 分布模糊集与最坏分布的联合求解技巧在Wasserstein DRO的双层max问题里最麻烦的是外层在“分布空间”上的搜索。理论上这个max可以写成对偶Wasserstein形式但有一个实践上很好用的转化利用对偶范数的性质把模糊集约束转化为对偶变量取值范围约束。简化地说当Wasserstein距离用1-范数或∞-范数定义时模糊集约束 \(W(Q, \hat{P}N) \le \epsilon\) 可以等价地写成关于辅助对偶变量 \(s_i\) 的约束 \[ \|s_i\|* \le \theta, \quad \sum_{i1}^N s_i \le \epsilon, \quad s_i \ge 0 \] 其中 \(\theta\) 是与真实风险参数相关的对偶标量\(s_i\) 可以理解为每个历史样本i的“惩罚价格”。这个转换的最大好处是原本在无穷维分布空间的max问题被压缩成了一组有限维的线性约束从而可以被常规优化求解器处理。用代码写这个过程的伪代码大概长这样# 伪代码Wasserstein DRO子问题的对偶转换 def build_sp(x_k, scenarios, eps, alpha): # scenarios: 历史样本 (w1, w2, ..., wN) # theta / s_i 是对偶变量 m OR.Model() s m.addVars(N, lb0) # 每个样本的惩罚对偶变量 theta m.addVar(lb0) # 风险对偶变量 # 模糊集约束 m.addConstr(quicksum(s[i] for i in range(N)) eps) for i in range(N): m.addConstr(s[i] -theta) # 对偶范数约束视范数类型而定 m.addConstr(s[i] theta) # 内层对偶目标 场景变量 z m.addVars(N, lb-GRB.INFINITY) # ... 这里拼接原问题的对偶目标 m.setObjective(quicksum(...) eps * theta, GRB.MAXIMIZE) return m上面只是一个轮廓实际写的时候要根据你的具体设备约束把目标函数和约束逐项写全。关键心得是不要试图直接求解“分布上的Max”一定要走对偶转化这条路径否则模型规模会膨胀到你根本算不动。3.4 算法收敛判据与初始可行解的选取CCG算法的迭代过程可以总结为两步循环主问题给出第一阶段决策 → 子问题找出最坏场景/分布并生成割平面 → 割平面加回主问题 → 再次求解主问题。这个循环一直持续到主问题的解不再显著变化。判断收敛通常用对偶间隙 \[ \text{Gap} \frac{\eta_{\text{MP}} - \text{SP}(x_k)}{|\text{SP}(x_k)|} \le \text{tol} \] 其中SP值为当前第一阶段决策下的最坏场景成本这是最优真实成本的上界主问题的目标值是当前算法的下界。工程上一般把 \(\text{tol}\) 设为0.1%或1%。设太紧会白白增加迭代轮数设太松又会得到一个不优的方案。实际测试下来0.1%到0.5%是比较合理的区间。初始可行解对迭代效率影响很大。为了避免第一轮主问题给出的方案完全不可行可以在主问题中加一些“启发式初始解”比如按照确定性预测值做一次最优机组组合把得到的启停状态作为初始可行解。这样第一轮子问题得到的割平面就不会太离谱整体迭代次数会明显减少。还有一个从实际调试里得来的经验每一轮割平面生成后不要把子问题的所有内部变量全保留在主问题里只保留当前轮的场景参数和对应调整变量就够。理论上每个割平面至少需要一组第二阶段变量但实际求解中场景数量太大如果都保留主问题的变量数量会疯狂膨胀到最后几十轮迭代之后MILP根本解不动。适当用“只保留每轮割平面相关变量”的策略可以在几乎不损失最优性的前提下大幅提升求解速度。4. 算例设计与对比分析4.1 测试系统基础参数与不确定场景生成为了验证模型效果我搭了一个简化但保留核心特征的独立微网测试系统火电机组1台最小出力30 MW最大出力150 MW爬坡率 ±3 MW/min最小启停时间4小时启停成本各5000元/次燃气轮机1台容量60 MW爬坡率 ±1.5 MW/min燃料成本比火电略高但碳排放低风电场装机80 MW光伏电站装机60 MW历史出力数据按30天、每15分钟一个点生成共2880个样本电池储能容量40 MWh最大功率20 MW充放电效率0.92/0.95SOC范围10%~90%电负荷预测曲线取典型工作日曲线峰值180 MW外网购电上限为80 MW实时购电价格约为日前价格的1.2倍。然后做三类对比模型确定性模型DET假设风电、光伏完全等于预测值没有不确定性传统box式鲁棒优化RO-BOX不确定参数在预测值正负20%范围内取最坏场景Wasserstein分布鲁棒优化DRO半径ε5%在经验分布周围建立模糊集DROCVaR置信度α0.95风险系数λ0.3在DRO基础上叠加尾部风险惩罚。这个对比其实就是一条“从乐观到悲观再到合理”的路线每个模型的差异完全体现在对不确定性的处理方式上。4.2 四类模型的横向对比结果我把四个模型跑完之后得到的结果很有意思也非常符合理论直觉模型日前成本(万元)实时期望调整成本(万元)总成本(万元)CVaR(万元)最坏场景成本(万元)DET85.218.6103.845.356.8RO-BOX102.45.2107.617.822.1DRO(ε5%)91.89.7101.527.434.2DROCVaR95.18.4103.518.926.5几个值得注意的结论第一确定性模型虽然昨日成本最低但把实时调整成本和最坏场景成本一算反而比所有鲁棒模型都要贵。这再次验证了在风光占比高的场景里忽略不确定性省下的钱最终都会在实时阶段连本带利赔回去。第二box式鲁棒模型的日前成本最高102.4万元因为它必须为最坏区间中的所有场景准备容量这是一个昂贵的保险。它的实时调整成本确实最低那只是表象——因为方案做得太保守系统空间被锁死了实际情况往往用不到那么多备用。第三DRO虽然日前成本比RO-BOX低10万以上实时阶段也没有过度恐慌总成本是最低的。原因很简单它不会像一个没见过世面的box集合那样把所有可能性都当成等概率的最坏情况而是根据历史样本分布特性做区分——那些离经验分布远到不合理的场景在模糊集那里就被过滤掉了。第四加CVaR之后日前成本上涨了3.3万元但CVaR从27.4万降到了18.9万下降了接近31%。这说明CVaR本质上是在“用低成本买高保障”对调度员来说非常划算。4.3 模糊集半径与风险系数的灵敏度分析半径ε的取值对最终结果的影响非常显著我实测的结果印证了理论推导。当ε从0开始增大时模型先从随机规划形态出发逐渐变得保守日前成本稳步上升实时期望成本逐渐下降。但到某个临界点之后上升和下降的斜率都趋于平缓——此时模型已经包含了几乎所有可能的分布再增大ε只是让日前成本白白上涨对安全性的改善微乎其微。实际调试中我习惯用第二节提到的 \(\epsilon_N(\beta)\) 公式先算出理论值然后在理论值附近做0.5倍、1倍、2倍三档灵敏度测试。如果运行结果在1倍和2倍之间差异不大说明模型对ε不敏感可以放心取1倍左右如果差异很大说明你的样本对不确定性的刻画能力不足可能需要增加历史数据而不是盲目加大半径。风险系数 \(\lambda\) 的规律也类似。\(\lambda\)从0增大到1.0的过程中CVaR会快速下降然后趋于饱和。我们在工程中关注的其实不是那个让CVaR降到最低的λ而是找到一个“性价比”最高点——再往上加风险厌恶系数成本上去了尾部风险却不怎么降了。这个拐点通常可以通过画出\lambda-CVaR曲线来寻找我自己在多个案例里得到的经验是最优拐点一般出现在0.2到0.5之间。5. 工程实现中的常见问题与避坑指南5.1 建模阶段的典型错误第一类错误也是最多人踩的坑忘记给子问题加合适的无界防护。我在复现之初子问题对偶之后经常出现unbounded排错排了整整两天最后发现是内层对偶中少了一个至关重要的松弛变量约束。功率平衡约束在数学上是一个等式约束如果不给对偶加入合适的边界比如允许弃风和切负荷的非负松弛变量那么当场景极端时对偶问题确实会出现无界。建议每一类设备都显式写入松弛变量哪怕你不想用它也要给它一个足够高的惩罚系数这样既保证了对偶有限性也符合实际调度的物理语义。第二类错误SOC递推约束跨阶段写串了。日前阶段的SOC会影响实时阶段的可用容量我在初版模型里把两阶段的SOC都写成一整条约束链结果实时阶段的最优方案把储能“从未来借电”来用算出来一个看似漂亮但不物理的结果。正确的做法是日前阶段决定SOC基准轨迹实时阶段在基准轨迹周围允许上下调整但必须在SOC上下限内重新约束。第三类错误风机和光伏取值超出物理可行域。Wasserstein模糊集半径设得太大时最坏分布的场景里可能出现负的风电出力或者超过装机容量的光伏出力。解决方法是在所有场景约束上显式加上物理上下界比如 \(0 \le \tilde{w}{wind} \le P{wind}^{cap}\)。这也是DRO模型中一个容易被忽略但非常必要的工程细节。5.2 求解阶段的数值问题数值问题是实际调bug时最磨人的部分。首先要注意量纲统一。我用过一套代码成本以万元为单位功率以MW为单位但储能SOC约束里出现千瓦和MWh混用结果目标函数波动幅度大得离谱gurobi直接报数值问题。建议从一开始就把所有物理量统一到MW和MWh成本单位用万元并保持全模型一致。其次是求解器的数值容差设置。对于这类包含大系数比如弃风惩罚5000元/MWh和普通系数燃料成本约300元/MWh的模型如果默认的BarHomogeneous、NumericFocus设置不够好会频繁出现数值警告。我在跑大规模场景时习惯把Gurobi的NumericFocus调到2同时把MIPGap设成0.5%这样最终结果的稳定性和可复现性都更好。子问题的Max-Max特性还会引入一个典型现象求出来的最坏场景往往位于模糊集或场景集合的边界上。这是正常的但如果你的求解器在每次迭代中给出的最坏场景在“边界上来回跳”说明割平面数目可能不够或者收敛容差设置得太松。曾经遇到过连续15次迭代目标值都没有变化但实际最优解差得远的怪事后来发现是子问题场景规模太大了求解器被MIPGap过早终止导致生成的割平面质量很差。解决办法是把子问题单独拿出来做一次精确求解检查它返回的SP值是否与承诺的收敛值一致。5.3 调参与收敛性经验迭代轮数方面在我的测试系统上CCG一般在15~30轮之间收敛到0.5%的间隙。如果你的模型跑到50轮以上还不收敛我建议从以下三个方向排查子问题的MIPGap是否太大导致割平面失真场景数量是否过多每轮加进主问题的变量规模超过了解题能力模糊集半径是否过大导致子问题的可行域异常复杂。还有一个很实用的技巧在迭代初期使用较大容差比如1%快速逼近可行域后面再收紧到0.1%做精细迭代。这样的“粗-精两阶段收敛”可以把总计算时间缩短约30%而且最终结果的精度几乎不受影响。5.4 发到CSDN这类技术社区时的一点心得体会虽然这篇博客的定位是技术向干货但如果你打算把模型复现的结果整理成文章发布有几个细节值得留意。一是代码结构要清晰。我推荐把“参数输入data.py”“模型构建model.py”“CCG迭代ccg_solver.py”“结果可视化plot.py”四个模块分开每份代码控制在200行以内。这样不仅方便读者复现也方便你自己后续的扩展。二是跑结果的截图要规范。这类优化调度模型的读者群体很多是研究生和工程师他们会非常关注目标函数值、迭代曲线和对比表格。建议保存好每个模型的Gap收敛曲线并且在文中明确标注求解器版本、硬件配置和求解时间——这些信息比“我们的方法很好”有说服力得多。三是把关键参数和结果以表格形式展示。无论你是放Github仓库链接还是在CSDN里贴完整代码一份清晰的结果对比表格能让文章在首页的曝光效率高很多。我个人在实际操作中的体会是Min-Max-Max-Min这种四层架构真正强大的地方不在于它比传统方法“多了一层Max”而在于它把“数据不确定性”和“场景不确定性”放进了同一个统一的优化框架里处理。工程上你可以用历史样本直接驱动模型不需要对风光分布做任何主观假设理论上它又能退化成随机规划和box鲁棒这两个极端给对比实验提供了很好的基准。最后再分享一个小技巧如果你手头的历史数据量不够大不要急着扩大Wasserstein半径来“保险”优先考虑把历史数据按天气类型分桶分别构造模糊集这样在总样本量不变的情况下每个桶内的分布特征会更清晰模型的保守度也有望更低。
