做微网调度优化的人多多少少都被不确定性搞过头疼。光伏和风电出力上一秒还在预测曲线上下一秒就可能断崖式下跌如果调度模型事先不知道会发生这种事那常规的确定性优化结果在实际运行时往往就变成废纸。我这次做的项目解决的正是这类问题基于关键场景辨别算法的两阶段鲁棒微网优化调度用Matlab从模型搭建到CCG迭代求解完整落地。项目思路不复杂——第一阶段做预调度第二阶段在不确定场景发生时做再调整而关键场景辨别算法负责从海量随机场景里挑出那些真正会让系统吃紧的“坏场景”避免枚举全部场景。整套代码我在Matlab R2023b YALMIP CPLEX环境下跑通了微网模型、场景生成与削减、主问题子问题迭代都能稳定工作。这篇博文我会把建模思路、算法细节、Matlab实现步骤和调试过程都拆开讲正在研究微网优化调度、主动配电网、鲁棒优化或者被CCG不收敛和场景太多算不动折磨的同学应该能从里面找到能直接用的方案。1. 项目读懂两阶段鲁棒微网优化调度到底在解决什么问题1.1 为什么确定性调度搞不定微网前两年我刚接触微网调度时习惯性用确定性优化预测明天光伏出力、风电出力、负荷曲线然后把所有预测值当成已知参数丢进优化模型解出一份“最优”的机组启停和储能充放电计划。听起来没什么问题但实际一跑就露馅——光伏在下午可能因为一片云就跌落一半负荷高峰也可能比预测提前两小时到来储能如果在上午按预测方案放光了电晚上面对突发高峰只能干瞪眼最后要么高价从主网购电要么直接切负荷。这种方案失效的本质原因是它把不确定性当成了确定性来处理。微网和大型电网不一样系统惯性小、可调资源有限、PCC联络线功率也有限预测误差一旦撞上运行约束的边界整个调度方案就会从“最优”变成“不可行”。所以要真正解决微网调度问题必须在模型里显式表达“未来可能发生什么”。这也是鲁棒优化的出发点接受预测会错并且为不确定参数建立集合让方案在最坏情况下也不要太难看。1.2 两阶段鲁棒把“后悔权”给了调度员单阶段鲁棒优化的想法很简单把所有决策都在不确定性实现之前做让整个方案在任何可能场景下都可行。这样做的结果往往非常保守因为微网里有大量“可以等看到实际天气再调整”的决策比如柴油机出力、储能充放电功率没必要在第一天就锁死。两阶段鲁棒优化把决策分成了两层。第一阶段决策是“现在就要拍板”的慢变量例如机组启停状态、与大电网的日前交互计划这类决策往往需要提前安排短时间内改不了。第二阶段决策是“看到实际场景后再说”的快变量例如在上述启停方案下各台机组多出力还是少出力、储能是充电还是放电、是否需要紧急切负荷。目标函数也相应变成第一阶段成本加上在所有可能场景中第二阶段调整成本最大的那一档。换句话说两阶段鲁棒给了调度员一个“后悔权”但同时要求他为最坏情况后的调整买单。这里可以打个比方提前订酒店和机票属于第一阶段到了当地发现下雨再决定是打车还是坐地铁属于第二阶段。模型不要求在出发前就把第二阶段的每一步想死而是允许在“天气实现”后再做调整这比单阶段鲁棒要灵活不少。1.3 关键场景辨别算法在这一套里扮演什么角色两阶段鲁棒模型的难点是第二阶段里那个 max-min 结构外层要在不确定集里找最恶劣的不确定参数内层要在该参数下找最小调整成本。理论上不确定集是连续空间组合无限实际工程里通常用采样生成大量离散场景来近似这个连续空间。但问题来了——如果每个CCG迭代都要对几千甚至几万个场景做精确子问题求解计算时间根本扛不住。关键场景辨别算法就是解决这个痛点的。它做的事情不是把场景一视同仁地削减而是通过某种评价指标优先识别出那些“可能让系统付出最高代价”的场景再让CCG只针对这些关键场景构造割约束。这样既保住了最坏情况又极大压缩了计算规模。可以把它理解为两阶段鲁棒框架里的加速器先从一堆随机样本里找出最有可能“砸场子”的选手再把训练重点放在对付这些选手上。我这次实现的基本流程是先生成大量候选场景做一轮快速场景削减然后在CCG迭代过程中用当前第一阶段解对削减后场景做快速评估筛选出Top K个关键场景进入精确子问题。整体计算下来场景规模从几万压缩到几十单次迭代耗时从分钟级降到秒级。2. 模型怎么搭不确定集、目标函数和约束条件逐个说清2.1 先确定你有哪些决策权和资源写代码之前得先把优化模型的边界画清楚。我这次用的是典型孤岛可并网微网结构包含柴油发电机、风电、光伏、储能电池以及不可控负荷通过PCC与上级电网相连既能购电也能售电。调度周期取24小时时间分辨率为1小时。这个粒度已经能覆盖大多数研究场景不用太细否则计算量会翻倍。决策变量怎么分层是我踩了很多坑后才理清的。第一阶段变量我放的是柴油机启停状态0-1变量、并网/孤岛运行状态、与大电网的日前购售电计划。第二阶段变量放的是各台柴油机的实际出力调整量、储能充放电功率、弃风弃光量、切负荷量。核心原则是0-1变量尽量留在第一阶段因为第二阶段如果有整数变量内层min问题就不是LP无法直接通过对偶转单层CCG实现会非常痛苦。如果你非要让储能充放电状态也在第二阶段那子问题就得处理成MILP事情复杂太多调试周期也会指数级上升。2.2 不确定集用盒式还是预算约束不确定参数我选了三个风电出力、光伏出力、节点负荷。它们都是微网里最容易预测失准的变量。不确定集最简单的是盒式集合也就是每个时刻的不确定参数在预测值上下一个固定偏差区间内浮动。但纯盒式太极端因为它允许所有不确定参数同时取最大偏差现实中哪有那么巧的事这种集合出来的结果往往保守得离谱。所以我用了带预算约束的不确定集。以风电为例设预测出力为 (P_{w,base}(t))最大预测偏差为 (d_w(t))实际出力可以写成(P_w(t) P_{w,base}(t) \delta_w(t) \cdot d_w(t))其中 (\delta_w(t)) 是一个0到1之间的不确定变量表示t时刻实际出力偏离预测的程度。为了控制总体的偏离程度再加一个预算约束(\sum_{t1}^{24} \delta_w(t) \leq \Gamma_w)(\Gamma_w) 就是预算值它允许在整个调度周期内最多只有若干时段出现较大偏差。这个参数非常关键(\Gamma_w) 越大模型越保守方案越贵但越抗风险(\Gamma_w) 设为0就退化成确定性模型。调试时我喜欢先设一个适中的 (\Gamma)比如12也就是允许一半时段出现偏差观察结果合理性再调。2.3 目标函数与约束条件的建模要点目标函数我分两部分。第一部分是第一阶段成本表达式如下(C_1 \sum_{t} \left[ \sum_{g} (c_{g}^{fuel} \cdot P_{g}(t) c_{g}^{su} \cdot u_{g}^{su}(t)) c_{buy}(t) \cdot P_{buy}(t) - c_{sell}(t) \cdot P_{sell}(t) \right])其中 (c_{g}^{fuel}) 是机组燃料成本系数(c_{g}^{su}) 是启动费用(u_{g}^{su}(t)) 是启动动作标识(c_{buy}/c_{sell}) 是购售电价(P_{buy}/P_{sell}) 是购售电功率。第二部分是第二阶段最坏情况下的调整成本包含柴油机额外出力费用、储能磨损成本、弃风弃光惩罚和切负荷惩罚。注意切负荷惩罚系数一定要设得比其他成本高一个数量级否则优化器会为了省钱而偷偷切负荷。约束条件我分成几组写代码时也是按组添加不然排查问题会很乱功率平衡约束任意时刻所有机组出力加储能放电加购电功率必须等于负荷加储能充电加售电功率。这是最基本也是最容易漏的约束。柴油机约束出力上下限、爬坡约束、最小启停时间约束。爬坡约束在第二阶段非常重要因为快速调整能力直接决定了鲁棒性。储能约束SOC递推式、充放电功率上下限、不能同时充放电。不能同时充放电可以用两个0-1变量加big-M处理但如果把第二阶段充放电状态也设成整数就破坏了LP结构。我这次的处理是把充放电状态放到第一阶段第二阶段只留功率连续变量。联络线约束购售电功率不能同时出现且不超过最大交换功率。不确定集约束如2.2节所述每个不确定变量在偏差范围内且满足预算约束。还有一个容易忽略的点第二阶段里柴油机出力调整是有物理边界的不能超过第一阶段机组状态对应的爬坡空间。这个约束我在初版代码里漏了导致子问题求出的调整成本偏低CCG上下界差距始终下不去查了很久才找到。2.4 第二阶段子问题为什么要转成单层两阶段鲁棒的子问题结构是给定第一阶段解 (x)求(\max_{u \in U} \min_{y \in F(x,u)} c^T y)也就是先在不确定集里找最恶劣的场景 (u)再在该场景下求最小调整成本。这个 max-min 结构看起来只有两层但在YALMIP里没法直接嵌套求解因为 (U) 也是优化变量内层min又是个优化问题整体是一个双层规划。标准做法是把内层min问题线性规划化并写出对偶问题。由于第二阶段变量全是连续的内层min是LP强对偶成立可以转换为一个max问题。然后这个max问题再和外层max合并变成单层最大化的线性规划问题。转换后决策变量包括不确定变量 (u) 和内层对偶变量 (\lambda)目标函数变成类似 ((b - B u)^T \lambda) 的形式。这个时候再交给YALMIP就只是一个带线性约束的最大化问题CPLEX轻松解出。我在调试时反复检查过转换过程里最容易错的是对偶约束符号、目标函数中 (u) 与 (\lambda) 的耦合项以及可行域 (F(x,u)) 里哪些约束需要保留到对偶中。这些细节一旦错一个算出来的“最恶劣场景”就是错的整个CCG迭代也会跟着跑偏。3. 核心算法拆解关键场景辨别 CCG迭代3.1 场景生成别用蒙特卡洛硬扛做随机规划的人第一反应就是蒙特卡洛随机采样简单是简单但有个天然问题为了覆盖小概率极端情况需要抽上万个样本否则最恶劣场景很可能压根没被抽到。而样本越多后续排序和精确子问题求解就越慢。我的建议是上采样阶段用拉丁超立方采样LHS至少比纯随机均匀得多。在Matlab里用lhsdesign(N, K)先生成N行K列的[0,1]均匀分布归一化样本再用逆累计分布函数映射到目标分布。风电和光伏出力我用Beta分布负荷用正态分布。比如光伏出力的采样可以用betainv(X, alpha, beta)乘以装机容量得到。这里的 alpha 和 beta 要根据实际预测偏差的历史数据去做最大似然估计实在没有历史数据就先设成让方差和预测偏差匹配的经验值。我实际的参数配置是初始采样5000个场景每个场景包含24小时的风电、光伏、负荷序列。注意这一步不能直接在内存里生成一个 5000×24×3 的密集矩阵然后开始枚举效率很低。更合理的是生成后立刻计算关键统计特征或者直接进入场景削减环节。3.2 场景削减把20000个场景浓缩成20个场景削减的目标是在尽量保留原始概率分布信息的前提下把场景数量从几千压到几十。常用的有K-means聚类和同步回代消除法。K-means实现简单但在处理场景概率时不够精细同步回代法更贴合随机规划的习惯它会维护每个场景的概率权重每次把两个距离最近的场景合并把其中一个的概率加到另一个上再更新距离矩阵直到达到目标场景数。我在Matlab里写了一个轻量版同步回代函数核心逻辑如下用pdist2计算所有场景两两之间的距离距离度量用欧氏距离不同变量之间先做归一化。找到距离最近的一对场景假设是场景i和场景j。删除其中一个场景同时把它的概率加到另一个场景上。更新时间点上的场景集合并重新计算距离。重复直到只剩目标数量。这个方法虽然理论推导看着复杂但工程实现很直接。唯一要注意的是场景按24小时序列比较时如果风电、光伏、负荷的量纲差别很大必须先将每个变量除以自己的基准值否则距离会被量纲大的变量主导削减出来一堆“伪相似”场景。把5000个场景削减到30个左右既能保留主要概率形态又能让后续关键场景辨别和CCG迭代跑得动。削减完一定要画图看看场景包络线如果削减后的场景集和原始预测曲线差太多那就是距离矩阵或者归一化出了问题。3.3 关键场景辨别寻找能让调度崩掉的场景到这里我要专门讲一下关键场景辨别算法这是整个项目里最容易被忽视但又最影响性能的部分。场景削减解决的是“数量多”的问题关键场景辨别解决的是“哪些场景才是真正危险的”的问题。什么意思呢两阶段鲁棒的CCG子问题本质上是在不确定集里找一个让第二阶段调整成本最高的 (u)。如果直接用连续不确定集求解子问题是个对偶后的LP但如果我们已经有一批离散场景也可以用一种更工程化的方式去逼近最恶劣场景先把每个候选场景固定下来分别代入当前第一阶段解 (x)计算该场景下第二阶段的最小运行成本然后按成本降序排序取前K个场景作为候选关键场景。为什么这样可以因为第二阶段成本越高说明这个场景对当前解的冲击越大。把这些高成本场景挑出来交给CCG去生成割约束效果上等同于在和系统最过不去的几个场景上做对抗。排序可以用近似模型先过滤比如少考虑几个约束只求目标值的一个上界然后在Top-N里再用完整子问题精确求解。这样能减少大量不必要的精确LP求解。我在代码里实现了一个key_scene_select.m它的输入是当前第一阶段解和候选场景集输出是排序后排名靠前的关键场景及其评价值。每次CCG迭代都会重新调用一次因为第一阶段解变了场景的危险程度排序也会变。这就是“动态关键场景辨别”比一开始就固定选几个场景要稳得多。3.4 CCG主从迭代的完整流程现在把CCG的迭代流程串起来。整个算法在项目里的执行顺序如下初始化设置最大迭代次数iter_max 20收敛阈值tol 1e-3下界LB -inf上界UB inf设置一个初始恶劣场景集合比如直接取预测场景。构建主问题主问题包含第一阶段决策变量以及针对当前已发现场景集合中的每个场景对应一组第二阶段变量和约束。目标函数是第一阶段成本加一个辅助变量 eta表示当前已发现场景中最坏情况下的第二阶段成本。求解主问题得到第一阶段最优解 (x_k) 和当前下界 (LB 第一阶段成本 eta)。固定 (x_k)调用关键场景辨别先对候选场景排序筛选Top K个候选场景然后对这些候选场景用精确子问题求最恶劣场景 (u^*)并得到子问题目标值 (f_k)。更新上界(UB \min(UB, 第一阶段成本(x_k) f_k))。如果 (UB - LB) 小于阈值跳出循环。否则把找出的最恶劣场景 (u^*) 加入主问题的场景集合并在主问题中新增一组第二阶段变量和约束也就是新增一列回到第3步。这里要特别注意主问题的规模会随着迭代次数增加而线性增长因为每个新增场景都带一套第二阶段变量。如果迭代20次主问题里就有20组第二阶段约束规模还行但如果迭代50次还没收敛要么约束写错要么场景辨别逻辑有问题主问题本身会越来越难解最终可能自己先把内存耗尽。4. Matlab代码实现从主程序到子问题一步步落地4.1 代码文件结构和运行前准备我这次项目的代码目录如下结构比较清晰可以供参考main_robust_schedule.m主函数负责参数初始化、调用场景处理和CCG循环。load_data.m载入微网参数比如柴油机爬坡率、储能容量、负荷曲线、风光预测曲线。scenario_generate.m生成初始场景集使用LHS采样。scenario_reduce.m同步回代法削减场景。key_scene_select.m关键场景辨别对候选场景排序筛选。build_master.m构建主问题的YALMIP模型并求解。build_subproblem.m构建给定第一阶段解下的精确子问题返回最恶劣场景和子问题目标值。运行环境方面Matlab版本建议R2020以上我用的是R2023b。YALMIP是必须的求解器推荐CPLEX或Gurobi两者都可以通过学术免费申请。如果暂时没有商业求解器也可以用intlinprog硬跑但性能会差很多尤其是主问题规模一大基本等不起。4.2 主循环上下界迭代与终止条件CCG主循环是整段代码的骨架我在main_robust_schedule.m里写了一个比较典型的迭代框架代码结构如下% 初始化 iter 1; LB -inf; UB inf; tol 1e-3; iter_max 20; % 初始场景集先给一个预测场景 scene_set P_predict; while iter iter_max % 构建并求解主问题 [x_k, master_obj, eta] build_master(scene_set, params); LB master_obj; % 关键场景辨别先生成大量候选场景削减再排序 cand_scenes scenario_generate(params.N_samples, params); cand_scenes scenario_reduce(cand_scenes, params.N_reduce); [key_scenes, sub_obj] key_scene_select(x_k, cand_scenes, params, params.TopK); % 用精确子问题求最恶劣场景 [u_star, f_k] build_subproblem(x_k, params); % 更新上界 UB min(UB, first_stage_cost(x_k, params) f_k); fprintf(iter %d, LB %.4f, UB %.4f, gap %.4f\n, ... iter, LB, UB, (UB - LB)/UB); if (UB - LB) / abs(UB) tol break; end % 把最恶劣场景加入主问题场景集合 scene_set [scene_set, u_star]; iter iter 1; % 这里还可以加一个限制如果scene_set容量超过预设最大值停止 end这段代码看着简单但有几个细节需要特别注意。第一build_master每次都要重新构建整个YALMIP模型因为我看到有人习惯把模型定义在循环外面、只更新约束结果经常出现变量维度不匹配。第二scene_set作为元胞数组或者纵向拼接的大矩阵新增场景时要注意维度我习惯把它定义成一个三维数组每一页是一个场景方便主问题里用循环展开。key_scene_select返回的sub_obj是候选场景的评估成本主要用来更新排序不一定直接参与CCG上界更新。真正用于上界更新的是build_subproblem得到的f_k因为那是通过精确对偶子问题求出来的最恶劣场景成本。初版代码里我把这两个搞混了导致UB和LB对不上迭代乱跳后来打印每一步的值才定位到问题。4.3 YALMIP建模与求解器配置主问题和子问题的建模我都是用YALMIP因为它的表达式写法接近数学公式调试起来直观很多。主问题的变量定义大致如下% 第一阶段变量 u_start binvar(n_gen, T, full); % 机组启停状态 u_su binvar(n_gen, T, full); % 启动动作 P_buy sdpvar(1, T, full); % 购电功率 P_sell sdpvar(1, T, full); % 售电功率 % 辅助变量 eta eta sdpvar(1, 1); % 对当前场景集中的每个场景定义一组第二阶段变量 for k 1:size(scene_set, 3) P_gen{k} sdpvar(n_gen, T, full); % 该场景下机组出力 P_ch{k} sdpvar(1, T, full); % 储能充电 P_dis{k} sdpvar(1, T, full); % 储能放电 P_curtail{k} sdpvar(2, T, full); % 弃风弃光 P_loadcut{k} sdpvar(1, T, full); % 切负荷 end这里有一个容易犯的错误第二阶段变量是在循环里随着新增场景不断增加的因此变量定义的维度要动态变化。YALMIP支持元胞数组存变量但如果你反复对同一个元胞数组扩展性能会下降。我习惯每次进入build_master时根据scene_set的当前大小重新生成全部变量这样代码清晰也不怕维度对不上。求解器配置方面我在主程序开头就设置了通用配置options sdpsettings(solver, cplex, verbose, 1, ... showprogress, 1, cplex.mip.tolerances.mipgap, 1e-4);如果机器上装的是Gurobi只需要把solver改成gurobi。注意YALMIP对求解器的名字很敏感大小写也要匹配不确定时可以用yalmiptest查看可用求解器列表。我在早期调试时经常因为求解器识别不了而报错后来都是先跑一次yalmiptest确认环境。4.4 子问题与制约生成的代码要点子问题的实现是整个项目里最核心也最容易出错的环节。我前面说过子问题是 max-min 结构首先要对第二阶段内层LP取对偶把问题变成单层max。这一步如果用纯代码实现最简单的做法是把内层约束写成矩阵形式 (A y \leq b B u)然后利用YALMIP定义对偶变量和对偶约束。内层问题的形式可以写成(\min_y \ c^T y)s.t. (A y \leq b B u)对偶变量设为 (\lambda \geq 0)对偶问题为(\max_{\lambda \geq 0} \ (b B u)^T \lambda)s.t. (A^T \lambda c)然后外层再对 (u) 最大化。合并后子问题的决策变量是 (u) 和 (\lambda)约束是 (u \in U)(\lambda \geq 0)(A^T \lambda c)目标函数是 ((b B u)^T \lambda)。在YALMIP里写出来大概是% 子问题变量不确定变量 u对偶变量 lambda u sdpvar(dim_u, 1, full); lambda sdpvar(dim_lambda, 1, full); % 对偶可行约束 constraints [A * lambda c, lambda 0]; % 不确定集约束 constraints [constraints, u_min u u_max]; constraints [constraints, sum(u_dev) budget]; % 目标函数 objective (b B * u) * lambda; % 求解最大化子问题 optimize(constraints, -objective, options); u_star value(u); f_k value(objective);这段代码看起来逻辑通顺但实际跑起来很容易因为A、B、b、c的维度对应不上而出错。我自己的经验是先写一个固定场景的内层LP用YALMIP求一次min再用对偶变量的dual方法检查对偶目标是否等于原目标。确认无误后再把 (u) 从固定值改成sdpvar用上面的合并形式求解。这样分两步走定位维度错误的效率高很多。另外要提一下可行性问题。如果在某个不确定参数下第二阶段问题不可行比如柴油机爬坡不够、储能容量不足那么对偶问题会无界或不可行。处理方法是在内层加入人工变量比如允许轻微切负荷、弃电并在目标函数里加上惩罚系数。惩罚系数设置要大于任何正常调整成本否则优化器会把切负荷当成一种常规手段而不是紧急手段。5. 调试与避坑我实际跑代码时踩过的坑5.1 求解器装不上或连不上怎么办很多第一次跑YALMIP的同学会遇到Solver not found或者No appropriate solver的报错。这个问题的根源不是YALMIP没识别到求解器就是Matlab路径里压根没有添加CPLEX或Gurobi的接口。CPLEX安装后要在Matlab里执行addpath(genpath(C:\Program Files\IBM\ILOG\CPLEX_Studio2211\cplex\matlab))然后保存路径。Gurobi类似要添加gurobi的Matlab接口目录。还有一个冷门坑CPLEX启动需要许可证服务有时候Matlab命令行里yalmiptest显示可用但一求解就卡住或报许可证错误。这种现象一般是许可证过期或者网络许可证连接不上重新启动许可证服务即可。如果实在解决不了商业求解器的问题小规模测试阶段可以用intlinprog应急。YALMIP支持把求解器设为intlinprog在optimize里会调用Matlab自带的整数线性规划求解器。缺点是速度慢但对于一个只有几台机组、24时段、十来个场景的小模型跑通流程还是没问题的。5.2 上下界长期不收敛的排查清单我在调试这个项目时最崩溃的阶段就是LB和UB死活不收敛。LB持续上涨UB却纹丝不动或者gap在某个值附近来回震荡。后来我总结了一套排查流程按顺序执行基本能定位问题先检查子问题目标值是否正常。固定第一阶段解随便挑一个极端场景手工算一下第二阶段成本和子问题输出对比。如果差别大说明子问题目标函数或对偶约束写错了。检查新加入主问题的场景是否真的是最恶劣场景。打印出 (u^*)看它是否在不确定集边界上。如果不是说明子问题搜索范围有问题可能是 (u) 被错误地限制在了一个小子集里。检查主问题中针对每个新场景添加的约束是否真正约束了该场景下的第二阶段变量。最容易犯的错误是复制上一段代码时约束里忘了把场景索引k换掉导致所有场景共用一套变量主问题变成了一个没有意义的约束LB自然上不去。如果gap一直在很小范围内震荡可能是收敛阈值设得太严格。两阶段鲁棒模型对数值误差比较敏感tol从1e-6放宽到1e-3往往就稳定了。另外不要一味依赖自动迭代。我在调试时会强制在第一次迭代后手动把预测场景和几个极端场景都加入主问题看主问题目标值是否明显上升。如果上升不明显说明这些场景并不是“恶劣场景”需要重新评估场景生成和场景辨别逻辑。5.3 性能优化和数值问题两阶段鲁棒最容易跑着跑着内存爆炸尤其当主问题场景数量增加时。我的性能优化经验有三条。第一减少不必要的精确子问题求解。关键场景辨别里排序阶段用近似模型做粗筛只对Top K做精确LP求解这个思路我前面提过实践下来能让单次迭代时间降一个数量级。第二在CCG主问题里YALMIP约束矩阵最好一次性拼好避免在循环里不断[constraints, constraints_new]。第三所有变量量纲尽量接近。如果电压、功率、储备容量混着用数值尺度差太大CPLEX内部会出数值问题导致不可行或收敛慢。我统一用了kW和kWh并经常检查约束矩阵的条件数。还有一个很实际的问题Matlab在长循环里频繁调用YALMIPoptimize很慢。我有时候会选择把CCG迭代次数控制在30次以内如果超过30次还没收敛查bug比继续跑循环更划算。5.4 常见问题速查表我在项目调试中整理了一份速查表遇到类似问题可以直接对号入座现象可能原因排查/解决办法YALMIP报Solver not found求解器未添加路径或未安装用yalmiptest检查添加CPLEX/Gurobi路径主问题求解时间暴涨场景数量过多、整数变量过多场景削减减少第一阶段0-1变量子问题不可行第二阶段可行域冲突加入人工变量和惩罚项上下界gap不降子问题对偶写错或主问题新场景约束失效检查对偶约束和目标函数做好场景索引最恶劣场景总在预测值附近不确定集预算约束设置过紧增大Γ或放宽偏差范围结果里切负荷量总是很大惩罚系数设置不合理将切负荷惩罚系数调到正常成本的10倍以上数值病态导致求解器报Infeasible量纲不统一或大M过大统一量纲减小big-M值这张表不能解决所有问题但大部分常见坑都能覆盖。很多时候项目跑不出来不是算法理论不对而是某个细节约束在代码里写歪了逐项排查效率远高过重新推公式。6. 这块代码怎么改造成你自己的项目6.1 从微网到配电网扩展如果你不是做微网而是做配电网优化调度这套框架也能平移。区别主要在第二阶段约束里要加入潮流方程。配电网一般用DistFlow方程描述功率流动方程本身是二次的但可以通过线性化或二阶锥松弛转换为线性凸问题。两阶段鲁棒里如果子问题变成二阶锥规划对偶性质没有LP那么干净但依然可以用凸对偶或通过求解器直接处理。我建议先做线性化版本跑通后再升级。关键场景辨别算法在配电网里更加重要因为配电网的节点数和场景维度都比微网大得多如果不用场景筛选光枚举负荷和光伏场景就能让求解器直接卡死。6.2 从两阶段鲁棒到分布式鲁棒如果你觉得两阶段鲁棒太保守可以考虑分布式鲁棒优化。它是在两阶段鲁棒的基础上不直接假设不确定参数最坏情况而是假设不确定参数的真实概率分布属于一个模糊集比如以Wasserstein距离定义的球。目标变成最小化在最坏概率分布下的期望成本。这种模型下关键场景辨别算法依然有用。你可以用历史数据生成大量场景再用场景削减和关键场景筛选得到支持集结合对偶变换求解。本质上它比纯两阶段鲁棒复杂一个层级但项目的基础框架——CCG迭代、场景处理、YALMIP建模——都不会白做。6.3 一点个人经验做这个项目前前后后折腾了大概两周最大感受是两阶段鲁棒优化很容易陷入“理论推导看起来很顺代码一跑全是问题”的尴尬。我的建议是不要一上来就追求复杂的微网模型。先用一台柴油机、一台储能、一个风电、一个负荷24小时只有预测场景和最恶劣场景两个场景把CCG循环跑通。一旦这个小例子能收敛再逐步增加机组的台数、加入光伏、扩大场景集。每加一个部件都要单独验证边界条件而不是把所有东西一次性堆进去不然排查bug的成本会指数上升。还有一点关键场景辨别算法单独拿出来讲可能不算特别新奇的理论但它确实能解决实际工程里“场景太多算不动”的真问题。把场景削减和动态排序做好CCG迭代速度能提升一个数量级。在做毕设或者写论文时这部分可以作为项目亮点来写也可以当作加速技巧重点是它是可以落地、可以被复现的。 基于关键场景辨别算法的两阶段鲁棒微网调度与CCG求解
做微网调度优化的人多多少少都被不确定性搞过头疼。光伏和风电出力上一秒还在预测曲线上下一秒就可能断崖式下跌如果调度模型事先不知道会发生这种事那常规的确定性优化结果在实际运行时往往就变成废纸。我这次做的项目解决的正是这类问题基于关键场景辨别算法的两阶段鲁棒微网优化调度用Matlab从模型搭建到CCG迭代求解完整落地。项目思路不复杂——第一阶段做预调度第二阶段在不确定场景发生时做再调整而关键场景辨别算法负责从海量随机场景里挑出那些真正会让系统吃紧的“坏场景”避免枚举全部场景。整套代码我在Matlab R2023b YALMIP CPLEX环境下跑通了微网模型、场景生成与削减、主问题子问题迭代都能稳定工作。这篇博文我会把建模思路、算法细节、Matlab实现步骤和调试过程都拆开讲正在研究微网优化调度、主动配电网、鲁棒优化或者被CCG不收敛和场景太多算不动折磨的同学应该能从里面找到能直接用的方案。1. 项目读懂两阶段鲁棒微网优化调度到底在解决什么问题1.1 为什么确定性调度搞不定微网前两年我刚接触微网调度时习惯性用确定性优化预测明天光伏出力、风电出力、负荷曲线然后把所有预测值当成已知参数丢进优化模型解出一份“最优”的机组启停和储能充放电计划。听起来没什么问题但实际一跑就露馅——光伏在下午可能因为一片云就跌落一半负荷高峰也可能比预测提前两小时到来储能如果在上午按预测方案放光了电晚上面对突发高峰只能干瞪眼最后要么高价从主网购电要么直接切负荷。这种方案失效的本质原因是它把不确定性当成了确定性来处理。微网和大型电网不一样系统惯性小、可调资源有限、PCC联络线功率也有限预测误差一旦撞上运行约束的边界整个调度方案就会从“最优”变成“不可行”。所以要真正解决微网调度问题必须在模型里显式表达“未来可能发生什么”。这也是鲁棒优化的出发点接受预测会错并且为不确定参数建立集合让方案在最坏情况下也不要太难看。1.2 两阶段鲁棒把“后悔权”给了调度员单阶段鲁棒优化的想法很简单把所有决策都在不确定性实现之前做让整个方案在任何可能场景下都可行。这样做的结果往往非常保守因为微网里有大量“可以等看到实际天气再调整”的决策比如柴油机出力、储能充放电功率没必要在第一天就锁死。两阶段鲁棒优化把决策分成了两层。第一阶段决策是“现在就要拍板”的慢变量例如机组启停状态、与大电网的日前交互计划这类决策往往需要提前安排短时间内改不了。第二阶段决策是“看到实际场景后再说”的快变量例如在上述启停方案下各台机组多出力还是少出力、储能是充电还是放电、是否需要紧急切负荷。目标函数也相应变成第一阶段成本加上在所有可能场景中第二阶段调整成本最大的那一档。换句话说两阶段鲁棒给了调度员一个“后悔权”但同时要求他为最坏情况后的调整买单。这里可以打个比方提前订酒店和机票属于第一阶段到了当地发现下雨再决定是打车还是坐地铁属于第二阶段。模型不要求在出发前就把第二阶段的每一步想死而是允许在“天气实现”后再做调整这比单阶段鲁棒要灵活不少。1.3 关键场景辨别算法在这一套里扮演什么角色两阶段鲁棒模型的难点是第二阶段里那个 max-min 结构外层要在不确定集里找最恶劣的不确定参数内层要在该参数下找最小调整成本。理论上不确定集是连续空间组合无限实际工程里通常用采样生成大量离散场景来近似这个连续空间。但问题来了——如果每个CCG迭代都要对几千甚至几万个场景做精确子问题求解计算时间根本扛不住。关键场景辨别算法就是解决这个痛点的。它做的事情不是把场景一视同仁地削减而是通过某种评价指标优先识别出那些“可能让系统付出最高代价”的场景再让CCG只针对这些关键场景构造割约束。这样既保住了最坏情况又极大压缩了计算规模。可以把它理解为两阶段鲁棒框架里的加速器先从一堆随机样本里找出最有可能“砸场子”的选手再把训练重点放在对付这些选手上。我这次实现的基本流程是先生成大量候选场景做一轮快速场景削减然后在CCG迭代过程中用当前第一阶段解对削减后场景做快速评估筛选出Top K个关键场景进入精确子问题。整体计算下来场景规模从几万压缩到几十单次迭代耗时从分钟级降到秒级。2. 模型怎么搭不确定集、目标函数和约束条件逐个说清2.1 先确定你有哪些决策权和资源写代码之前得先把优化模型的边界画清楚。我这次用的是典型孤岛可并网微网结构包含柴油发电机、风电、光伏、储能电池以及不可控负荷通过PCC与上级电网相连既能购电也能售电。调度周期取24小时时间分辨率为1小时。这个粒度已经能覆盖大多数研究场景不用太细否则计算量会翻倍。决策变量怎么分层是我踩了很多坑后才理清的。第一阶段变量我放的是柴油机启停状态0-1变量、并网/孤岛运行状态、与大电网的日前购售电计划。第二阶段变量放的是各台柴油机的实际出力调整量、储能充放电功率、弃风弃光量、切负荷量。核心原则是0-1变量尽量留在第一阶段因为第二阶段如果有整数变量内层min问题就不是LP无法直接通过对偶转单层CCG实现会非常痛苦。如果你非要让储能充放电状态也在第二阶段那子问题就得处理成MILP事情复杂太多调试周期也会指数级上升。2.2 不确定集用盒式还是预算约束不确定参数我选了三个风电出力、光伏出力、节点负荷。它们都是微网里最容易预测失准的变量。不确定集最简单的是盒式集合也就是每个时刻的不确定参数在预测值上下一个固定偏差区间内浮动。但纯盒式太极端因为它允许所有不确定参数同时取最大偏差现实中哪有那么巧的事这种集合出来的结果往往保守得离谱。所以我用了带预算约束的不确定集。以风电为例设预测出力为 (P_{w,base}(t))最大预测偏差为 (d_w(t))实际出力可以写成(P_w(t) P_{w,base}(t) \delta_w(t) \cdot d_w(t))其中 (\delta_w(t)) 是一个0到1之间的不确定变量表示t时刻实际出力偏离预测的程度。为了控制总体的偏离程度再加一个预算约束(\sum_{t1}^{24} \delta_w(t) \leq \Gamma_w)(\Gamma_w) 就是预算值它允许在整个调度周期内最多只有若干时段出现较大偏差。这个参数非常关键(\Gamma_w) 越大模型越保守方案越贵但越抗风险(\Gamma_w) 设为0就退化成确定性模型。调试时我喜欢先设一个适中的 (\Gamma)比如12也就是允许一半时段出现偏差观察结果合理性再调。2.3 目标函数与约束条件的建模要点目标函数我分两部分。第一部分是第一阶段成本表达式如下(C_1 \sum_{t} \left[ \sum_{g} (c_{g}^{fuel} \cdot P_{g}(t) c_{g}^{su} \cdot u_{g}^{su}(t)) c_{buy}(t) \cdot P_{buy}(t) - c_{sell}(t) \cdot P_{sell}(t) \right])其中 (c_{g}^{fuel}) 是机组燃料成本系数(c_{g}^{su}) 是启动费用(u_{g}^{su}(t)) 是启动动作标识(c_{buy}/c_{sell}) 是购售电价(P_{buy}/P_{sell}) 是购售电功率。第二部分是第二阶段最坏情况下的调整成本包含柴油机额外出力费用、储能磨损成本、弃风弃光惩罚和切负荷惩罚。注意切负荷惩罚系数一定要设得比其他成本高一个数量级否则优化器会为了省钱而偷偷切负荷。约束条件我分成几组写代码时也是按组添加不然排查问题会很乱功率平衡约束任意时刻所有机组出力加储能放电加购电功率必须等于负荷加储能充电加售电功率。这是最基本也是最容易漏的约束。柴油机约束出力上下限、爬坡约束、最小启停时间约束。爬坡约束在第二阶段非常重要因为快速调整能力直接决定了鲁棒性。储能约束SOC递推式、充放电功率上下限、不能同时充放电。不能同时充放电可以用两个0-1变量加big-M处理但如果把第二阶段充放电状态也设成整数就破坏了LP结构。我这次的处理是把充放电状态放到第一阶段第二阶段只留功率连续变量。联络线约束购售电功率不能同时出现且不超过最大交换功率。不确定集约束如2.2节所述每个不确定变量在偏差范围内且满足预算约束。还有一个容易忽略的点第二阶段里柴油机出力调整是有物理边界的不能超过第一阶段机组状态对应的爬坡空间。这个约束我在初版代码里漏了导致子问题求出的调整成本偏低CCG上下界差距始终下不去查了很久才找到。2.4 第二阶段子问题为什么要转成单层两阶段鲁棒的子问题结构是给定第一阶段解 (x)求(\max_{u \in U} \min_{y \in F(x,u)} c^T y)也就是先在不确定集里找最恶劣的场景 (u)再在该场景下求最小调整成本。这个 max-min 结构看起来只有两层但在YALMIP里没法直接嵌套求解因为 (U) 也是优化变量内层min又是个优化问题整体是一个双层规划。标准做法是把内层min问题线性规划化并写出对偶问题。由于第二阶段变量全是连续的内层min是LP强对偶成立可以转换为一个max问题。然后这个max问题再和外层max合并变成单层最大化的线性规划问题。转换后决策变量包括不确定变量 (u) 和内层对偶变量 (\lambda)目标函数变成类似 ((b - B u)^T \lambda) 的形式。这个时候再交给YALMIP就只是一个带线性约束的最大化问题CPLEX轻松解出。我在调试时反复检查过转换过程里最容易错的是对偶约束符号、目标函数中 (u) 与 (\lambda) 的耦合项以及可行域 (F(x,u)) 里哪些约束需要保留到对偶中。这些细节一旦错一个算出来的“最恶劣场景”就是错的整个CCG迭代也会跟着跑偏。3. 核心算法拆解关键场景辨别 CCG迭代3.1 场景生成别用蒙特卡洛硬扛做随机规划的人第一反应就是蒙特卡洛随机采样简单是简单但有个天然问题为了覆盖小概率极端情况需要抽上万个样本否则最恶劣场景很可能压根没被抽到。而样本越多后续排序和精确子问题求解就越慢。我的建议是上采样阶段用拉丁超立方采样LHS至少比纯随机均匀得多。在Matlab里用lhsdesign(N, K)先生成N行K列的[0,1]均匀分布归一化样本再用逆累计分布函数映射到目标分布。风电和光伏出力我用Beta分布负荷用正态分布。比如光伏出力的采样可以用betainv(X, alpha, beta)乘以装机容量得到。这里的 alpha 和 beta 要根据实际预测偏差的历史数据去做最大似然估计实在没有历史数据就先设成让方差和预测偏差匹配的经验值。我实际的参数配置是初始采样5000个场景每个场景包含24小时的风电、光伏、负荷序列。注意这一步不能直接在内存里生成一个 5000×24×3 的密集矩阵然后开始枚举效率很低。更合理的是生成后立刻计算关键统计特征或者直接进入场景削减环节。3.2 场景削减把20000个场景浓缩成20个场景削减的目标是在尽量保留原始概率分布信息的前提下把场景数量从几千压到几十。常用的有K-means聚类和同步回代消除法。K-means实现简单但在处理场景概率时不够精细同步回代法更贴合随机规划的习惯它会维护每个场景的概率权重每次把两个距离最近的场景合并把其中一个的概率加到另一个上再更新距离矩阵直到达到目标场景数。我在Matlab里写了一个轻量版同步回代函数核心逻辑如下用pdist2计算所有场景两两之间的距离距离度量用欧氏距离不同变量之间先做归一化。找到距离最近的一对场景假设是场景i和场景j。删除其中一个场景同时把它的概率加到另一个场景上。更新时间点上的场景集合并重新计算距离。重复直到只剩目标数量。这个方法虽然理论推导看着复杂但工程实现很直接。唯一要注意的是场景按24小时序列比较时如果风电、光伏、负荷的量纲差别很大必须先将每个变量除以自己的基准值否则距离会被量纲大的变量主导削减出来一堆“伪相似”场景。把5000个场景削减到30个左右既能保留主要概率形态又能让后续关键场景辨别和CCG迭代跑得动。削减完一定要画图看看场景包络线如果削减后的场景集和原始预测曲线差太多那就是距离矩阵或者归一化出了问题。3.3 关键场景辨别寻找能让调度崩掉的场景到这里我要专门讲一下关键场景辨别算法这是整个项目里最容易被忽视但又最影响性能的部分。场景削减解决的是“数量多”的问题关键场景辨别解决的是“哪些场景才是真正危险的”的问题。什么意思呢两阶段鲁棒的CCG子问题本质上是在不确定集里找一个让第二阶段调整成本最高的 (u)。如果直接用连续不确定集求解子问题是个对偶后的LP但如果我们已经有一批离散场景也可以用一种更工程化的方式去逼近最恶劣场景先把每个候选场景固定下来分别代入当前第一阶段解 (x)计算该场景下第二阶段的最小运行成本然后按成本降序排序取前K个场景作为候选关键场景。为什么这样可以因为第二阶段成本越高说明这个场景对当前解的冲击越大。把这些高成本场景挑出来交给CCG去生成割约束效果上等同于在和系统最过不去的几个场景上做对抗。排序可以用近似模型先过滤比如少考虑几个约束只求目标值的一个上界然后在Top-N里再用完整子问题精确求解。这样能减少大量不必要的精确LP求解。我在代码里实现了一个key_scene_select.m它的输入是当前第一阶段解和候选场景集输出是排序后排名靠前的关键场景及其评价值。每次CCG迭代都会重新调用一次因为第一阶段解变了场景的危险程度排序也会变。这就是“动态关键场景辨别”比一开始就固定选几个场景要稳得多。3.4 CCG主从迭代的完整流程现在把CCG的迭代流程串起来。整个算法在项目里的执行顺序如下初始化设置最大迭代次数iter_max 20收敛阈值tol 1e-3下界LB -inf上界UB inf设置一个初始恶劣场景集合比如直接取预测场景。构建主问题主问题包含第一阶段决策变量以及针对当前已发现场景集合中的每个场景对应一组第二阶段变量和约束。目标函数是第一阶段成本加一个辅助变量 eta表示当前已发现场景中最坏情况下的第二阶段成本。求解主问题得到第一阶段最优解 (x_k) 和当前下界 (LB 第一阶段成本 eta)。固定 (x_k)调用关键场景辨别先对候选场景排序筛选Top K个候选场景然后对这些候选场景用精确子问题求最恶劣场景 (u^*)并得到子问题目标值 (f_k)。更新上界(UB \min(UB, 第一阶段成本(x_k) f_k))。如果 (UB - LB) 小于阈值跳出循环。否则把找出的最恶劣场景 (u^*) 加入主问题的场景集合并在主问题中新增一组第二阶段变量和约束也就是新增一列回到第3步。这里要特别注意主问题的规模会随着迭代次数增加而线性增长因为每个新增场景都带一套第二阶段变量。如果迭代20次主问题里就有20组第二阶段约束规模还行但如果迭代50次还没收敛要么约束写错要么场景辨别逻辑有问题主问题本身会越来越难解最终可能自己先把内存耗尽。4. Matlab代码实现从主程序到子问题一步步落地4.1 代码文件结构和运行前准备我这次项目的代码目录如下结构比较清晰可以供参考main_robust_schedule.m主函数负责参数初始化、调用场景处理和CCG循环。load_data.m载入微网参数比如柴油机爬坡率、储能容量、负荷曲线、风光预测曲线。scenario_generate.m生成初始场景集使用LHS采样。scenario_reduce.m同步回代法削减场景。key_scene_select.m关键场景辨别对候选场景排序筛选。build_master.m构建主问题的YALMIP模型并求解。build_subproblem.m构建给定第一阶段解下的精确子问题返回最恶劣场景和子问题目标值。运行环境方面Matlab版本建议R2020以上我用的是R2023b。YALMIP是必须的求解器推荐CPLEX或Gurobi两者都可以通过学术免费申请。如果暂时没有商业求解器也可以用intlinprog硬跑但性能会差很多尤其是主问题规模一大基本等不起。4.2 主循环上下界迭代与终止条件CCG主循环是整段代码的骨架我在main_robust_schedule.m里写了一个比较典型的迭代框架代码结构如下% 初始化 iter 1; LB -inf; UB inf; tol 1e-3; iter_max 20; % 初始场景集先给一个预测场景 scene_set P_predict; while iter iter_max % 构建并求解主问题 [x_k, master_obj, eta] build_master(scene_set, params); LB master_obj; % 关键场景辨别先生成大量候选场景削减再排序 cand_scenes scenario_generate(params.N_samples, params); cand_scenes scenario_reduce(cand_scenes, params.N_reduce); [key_scenes, sub_obj] key_scene_select(x_k, cand_scenes, params, params.TopK); % 用精确子问题求最恶劣场景 [u_star, f_k] build_subproblem(x_k, params); % 更新上界 UB min(UB, first_stage_cost(x_k, params) f_k); fprintf(iter %d, LB %.4f, UB %.4f, gap %.4f\n, ... iter, LB, UB, (UB - LB)/UB); if (UB - LB) / abs(UB) tol break; end % 把最恶劣场景加入主问题场景集合 scene_set [scene_set, u_star]; iter iter 1; % 这里还可以加一个限制如果scene_set容量超过预设最大值停止 end这段代码看着简单但有几个细节需要特别注意。第一build_master每次都要重新构建整个YALMIP模型因为我看到有人习惯把模型定义在循环外面、只更新约束结果经常出现变量维度不匹配。第二scene_set作为元胞数组或者纵向拼接的大矩阵新增场景时要注意维度我习惯把它定义成一个三维数组每一页是一个场景方便主问题里用循环展开。key_scene_select返回的sub_obj是候选场景的评估成本主要用来更新排序不一定直接参与CCG上界更新。真正用于上界更新的是build_subproblem得到的f_k因为那是通过精确对偶子问题求出来的最恶劣场景成本。初版代码里我把这两个搞混了导致UB和LB对不上迭代乱跳后来打印每一步的值才定位到问题。4.3 YALMIP建模与求解器配置主问题和子问题的建模我都是用YALMIP因为它的表达式写法接近数学公式调试起来直观很多。主问题的变量定义大致如下% 第一阶段变量 u_start binvar(n_gen, T, full); % 机组启停状态 u_su binvar(n_gen, T, full); % 启动动作 P_buy sdpvar(1, T, full); % 购电功率 P_sell sdpvar(1, T, full); % 售电功率 % 辅助变量 eta eta sdpvar(1, 1); % 对当前场景集中的每个场景定义一组第二阶段变量 for k 1:size(scene_set, 3) P_gen{k} sdpvar(n_gen, T, full); % 该场景下机组出力 P_ch{k} sdpvar(1, T, full); % 储能充电 P_dis{k} sdpvar(1, T, full); % 储能放电 P_curtail{k} sdpvar(2, T, full); % 弃风弃光 P_loadcut{k} sdpvar(1, T, full); % 切负荷 end这里有一个容易犯的错误第二阶段变量是在循环里随着新增场景不断增加的因此变量定义的维度要动态变化。YALMIP支持元胞数组存变量但如果你反复对同一个元胞数组扩展性能会下降。我习惯每次进入build_master时根据scene_set的当前大小重新生成全部变量这样代码清晰也不怕维度对不上。求解器配置方面我在主程序开头就设置了通用配置options sdpsettings(solver, cplex, verbose, 1, ... showprogress, 1, cplex.mip.tolerances.mipgap, 1e-4);如果机器上装的是Gurobi只需要把solver改成gurobi。注意YALMIP对求解器的名字很敏感大小写也要匹配不确定时可以用yalmiptest查看可用求解器列表。我在早期调试时经常因为求解器识别不了而报错后来都是先跑一次yalmiptest确认环境。4.4 子问题与制约生成的代码要点子问题的实现是整个项目里最核心也最容易出错的环节。我前面说过子问题是 max-min 结构首先要对第二阶段内层LP取对偶把问题变成单层max。这一步如果用纯代码实现最简单的做法是把内层约束写成矩阵形式 (A y \leq b B u)然后利用YALMIP定义对偶变量和对偶约束。内层问题的形式可以写成(\min_y \ c^T y)s.t. (A y \leq b B u)对偶变量设为 (\lambda \geq 0)对偶问题为(\max_{\lambda \geq 0} \ (b B u)^T \lambda)s.t. (A^T \lambda c)然后外层再对 (u) 最大化。合并后子问题的决策变量是 (u) 和 (\lambda)约束是 (u \in U)(\lambda \geq 0)(A^T \lambda c)目标函数是 ((b B u)^T \lambda)。在YALMIP里写出来大概是% 子问题变量不确定变量 u对偶变量 lambda u sdpvar(dim_u, 1, full); lambda sdpvar(dim_lambda, 1, full); % 对偶可行约束 constraints [A * lambda c, lambda 0]; % 不确定集约束 constraints [constraints, u_min u u_max]; constraints [constraints, sum(u_dev) budget]; % 目标函数 objective (b B * u) * lambda; % 求解最大化子问题 optimize(constraints, -objective, options); u_star value(u); f_k value(objective);这段代码看起来逻辑通顺但实际跑起来很容易因为A、B、b、c的维度对应不上而出错。我自己的经验是先写一个固定场景的内层LP用YALMIP求一次min再用对偶变量的dual方法检查对偶目标是否等于原目标。确认无误后再把 (u) 从固定值改成sdpvar用上面的合并形式求解。这样分两步走定位维度错误的效率高很多。另外要提一下可行性问题。如果在某个不确定参数下第二阶段问题不可行比如柴油机爬坡不够、储能容量不足那么对偶问题会无界或不可行。处理方法是在内层加入人工变量比如允许轻微切负荷、弃电并在目标函数里加上惩罚系数。惩罚系数设置要大于任何正常调整成本否则优化器会把切负荷当成一种常规手段而不是紧急手段。5. 调试与避坑我实际跑代码时踩过的坑5.1 求解器装不上或连不上怎么办很多第一次跑YALMIP的同学会遇到Solver not found或者No appropriate solver的报错。这个问题的根源不是YALMIP没识别到求解器就是Matlab路径里压根没有添加CPLEX或Gurobi的接口。CPLEX安装后要在Matlab里执行addpath(genpath(C:\Program Files\IBM\ILOG\CPLEX_Studio2211\cplex\matlab))然后保存路径。Gurobi类似要添加gurobi的Matlab接口目录。还有一个冷门坑CPLEX启动需要许可证服务有时候Matlab命令行里yalmiptest显示可用但一求解就卡住或报许可证错误。这种现象一般是许可证过期或者网络许可证连接不上重新启动许可证服务即可。如果实在解决不了商业求解器的问题小规模测试阶段可以用intlinprog应急。YALMIP支持把求解器设为intlinprog在optimize里会调用Matlab自带的整数线性规划求解器。缺点是速度慢但对于一个只有几台机组、24时段、十来个场景的小模型跑通流程还是没问题的。5.2 上下界长期不收敛的排查清单我在调试这个项目时最崩溃的阶段就是LB和UB死活不收敛。LB持续上涨UB却纹丝不动或者gap在某个值附近来回震荡。后来我总结了一套排查流程按顺序执行基本能定位问题先检查子问题目标值是否正常。固定第一阶段解随便挑一个极端场景手工算一下第二阶段成本和子问题输出对比。如果差别大说明子问题目标函数或对偶约束写错了。检查新加入主问题的场景是否真的是最恶劣场景。打印出 (u^*)看它是否在不确定集边界上。如果不是说明子问题搜索范围有问题可能是 (u) 被错误地限制在了一个小子集里。检查主问题中针对每个新场景添加的约束是否真正约束了该场景下的第二阶段变量。最容易犯的错误是复制上一段代码时约束里忘了把场景索引k换掉导致所有场景共用一套变量主问题变成了一个没有意义的约束LB自然上不去。如果gap一直在很小范围内震荡可能是收敛阈值设得太严格。两阶段鲁棒模型对数值误差比较敏感tol从1e-6放宽到1e-3往往就稳定了。另外不要一味依赖自动迭代。我在调试时会强制在第一次迭代后手动把预测场景和几个极端场景都加入主问题看主问题目标值是否明显上升。如果上升不明显说明这些场景并不是“恶劣场景”需要重新评估场景生成和场景辨别逻辑。5.3 性能优化和数值问题两阶段鲁棒最容易跑着跑着内存爆炸尤其当主问题场景数量增加时。我的性能优化经验有三条。第一减少不必要的精确子问题求解。关键场景辨别里排序阶段用近似模型做粗筛只对Top K做精确LP求解这个思路我前面提过实践下来能让单次迭代时间降一个数量级。第二在CCG主问题里YALMIP约束矩阵最好一次性拼好避免在循环里不断[constraints, constraints_new]。第三所有变量量纲尽量接近。如果电压、功率、储备容量混着用数值尺度差太大CPLEX内部会出数值问题导致不可行或收敛慢。我统一用了kW和kWh并经常检查约束矩阵的条件数。还有一个很实际的问题Matlab在长循环里频繁调用YALMIPoptimize很慢。我有时候会选择把CCG迭代次数控制在30次以内如果超过30次还没收敛查bug比继续跑循环更划算。5.4 常见问题速查表我在项目调试中整理了一份速查表遇到类似问题可以直接对号入座现象可能原因排查/解决办法YALMIP报Solver not found求解器未添加路径或未安装用yalmiptest检查添加CPLEX/Gurobi路径主问题求解时间暴涨场景数量过多、整数变量过多场景削减减少第一阶段0-1变量子问题不可行第二阶段可行域冲突加入人工变量和惩罚项上下界gap不降子问题对偶写错或主问题新场景约束失效检查对偶约束和目标函数做好场景索引最恶劣场景总在预测值附近不确定集预算约束设置过紧增大Γ或放宽偏差范围结果里切负荷量总是很大惩罚系数设置不合理将切负荷惩罚系数调到正常成本的10倍以上数值病态导致求解器报Infeasible量纲不统一或大M过大统一量纲减小big-M值这张表不能解决所有问题但大部分常见坑都能覆盖。很多时候项目跑不出来不是算法理论不对而是某个细节约束在代码里写歪了逐项排查效率远高过重新推公式。6. 这块代码怎么改造成你自己的项目6.1 从微网到配电网扩展如果你不是做微网而是做配电网优化调度这套框架也能平移。区别主要在第二阶段约束里要加入潮流方程。配电网一般用DistFlow方程描述功率流动方程本身是二次的但可以通过线性化或二阶锥松弛转换为线性凸问题。两阶段鲁棒里如果子问题变成二阶锥规划对偶性质没有LP那么干净但依然可以用凸对偶或通过求解器直接处理。我建议先做线性化版本跑通后再升级。关键场景辨别算法在配电网里更加重要因为配电网的节点数和场景维度都比微网大得多如果不用场景筛选光枚举负荷和光伏场景就能让求解器直接卡死。6.2 从两阶段鲁棒到分布式鲁棒如果你觉得两阶段鲁棒太保守可以考虑分布式鲁棒优化。它是在两阶段鲁棒的基础上不直接假设不确定参数最坏情况而是假设不确定参数的真实概率分布属于一个模糊集比如以Wasserstein距离定义的球。目标变成最小化在最坏概率分布下的期望成本。这种模型下关键场景辨别算法依然有用。你可以用历史数据生成大量场景再用场景削减和关键场景筛选得到支持集结合对偶变换求解。本质上它比纯两阶段鲁棒复杂一个层级但项目的基础框架——CCG迭代、场景处理、YALMIP建模——都不会白做。6.3 一点个人经验做这个项目前前后后折腾了大概两周最大感受是两阶段鲁棒优化很容易陷入“理论推导看起来很顺代码一跑全是问题”的尴尬。我的建议是不要一上来就追求复杂的微网模型。先用一台柴油机、一台储能、一个风电、一个负荷24小时只有预测场景和最恶劣场景两个场景把CCG循环跑通。一旦这个小例子能收敛再逐步增加机组的台数、加入光伏、扩大场景集。每加一个部件都要单独验证边界条件而不是把所有东西一次性堆进去不然排查bug的成本会指数上升。还有一点关键场景辨别算法单独拿出来讲可能不算特别新奇的理论但它确实能解决实际工程里“场景太多算不动”的真问题。把场景削减和动态排序做好CCG迭代速度能提升一个数量级。在做毕设或者写论文时这部分可以作为项目亮点来写也可以当作加速技巧重点是它是可以落地、可以被复现的。