两阶段鲁棒优化微网调度：关键场景辨别算法与Matlab实现

两阶段鲁棒优化在微网调度里这几年是真火但很多人一上手就卡在“不确定性集合怎么建”“场景怎么选”“迭代求解怎么收敛”这几个坎上。我自己当初用Matlab做这个课题的时候也是把文献翻了个底朝天代码一行行啃才把整套流程跑通。这篇就把我基于关键场景辨别算法的两阶段鲁棒微网优化调度完整思路、数学模型和Matlab代码实现细节整理出来给正要入坑或者被坑得不浅的朋友一份能直接参考的实操笔记。1. 核心思路为什么微网调度需要两阶段鲁棒1.1 微网调度到底难在哪微网调度本质上是一个经济调度问题目标就是在满足用户负荷需求的前提下让光伏、风电、储能、微型燃气轮机这些分布式电源协调运行使得总运行成本最低。但这里有个天然的麻烦光伏出力和负荷需求都是不确定的天有不测风云这句话用在微网调度上再合适不过。光伏今天能发100 kW明天可能一片云飘过来就只有30 kW负荷也是商业区和居民区的用电曲线差异很大同一个时段的负荷预测值跟实际值之间总存在偏差。传统确定性调度是把这些预测值当成真实值来算得到一个调度方案实际上执行的时候如果光伏突然掉一半方案就废了该切负荷的切负荷该买电的买电成本直接起飞。所以这几年大家把目光转向了鲁棒优化——我按最坏情况来做决策不管光伏怎么波动、负荷怎么变我的方案都不会导致系统崩溃或者成本失控。1.2 两阶段鲁棒的直观理解两阶段鲁棒优化用大白话说就是我今天先做一个“现在必须定下来”的决定比如机组开停机状态、和主网的购售电协议这些属于第一阶段决策也叫here-and-now决策。等明天真实的光伏出力和负荷数据出来了我再根据实际场景去调整“当时可以灵活变”的量比如储能充放电功率、微型燃气轮机的出力、从主网购电的功率这些属于第二阶段决策也叫wait-and-see决策。第二阶段的存在就是鲁棒性的核心。不确定性在第二阶段以最恶劣的形式出现我不管它怎么变第二阶段总有办法用可调变量去应对只要应对得住整个系统就是安全的。这个逻辑本质上就是先做好最坏情况下的预案然后等真实情况来了见招拆招。1.3 关键场景辨别算法解决什么问题标准的两阶段鲁棒优化是一个min-max-min三层结构直接求解非常困难。常见做法是CCG列与约束生成算法通过主子问题迭代把最恶劣场景找出来加入模型中。但CCG有个痛点每次迭代都要解一个第二阶段max-min子问题如果微网节点多、机组多、时间尺度细这个子问题的求解非常耗时迭代几十次才收敛是常有的事。关键场景辨别算法的思路就是在迭代开始之前或者迭代过程中先从不确定集合里筛选出少量“关键场景”——这些场景对目标函数影响最大几乎能代表整个不确定集合的最坏情况——只用这些场景参与迭代就能用很小的计算代价换来接近全集合鲁棒的解。可以把它理解为“先侦察敌情再集中火力打关键目标”而不是把整个战场无差别轰炸一遍。注意关键场景的筛选是有讲究的不是简单随机抽几个场景。选少了解不够鲁棒选多了计算量又上去了。实际中常用的方法是结合对偶变量信息或者场景聚类来做这个后面在数学模型和代码部分会详细拆。2. 数学模型目标函数、约束与不确定性集合详解2.1 第一阶段目标函数与决策变量我们考虑一个典型的交流微网包含光伏PV、风电WT、微型燃气轮机MT、储能ESS以及与主网的联络线。调度的时域取24小时时间间隔1小时这样调度模型规模适中也是论文里最常见的配置。第一阶段决策变量主要是机组启停状态以及是否与主网签订购售电协议。目标函数的第一阶段部分包括微型燃气轮机的启动成本每次开机都要付一笔启动费类似于你打车起步价。如果考虑了开停机状态相关的固定运行成本也放在这一阶段。第一阶段目标可以写成[ \min_{\mathbf{x}} \left( \sum_{t} \sum_{g} SU_{g} \cdot y_{g,t} C_{fixed}(\mathbf{x}) \max_{\mathbf{u} \in \mathcal{U}} \min_{\mathbf{y} \in \mathcal{F}(\mathbf{x}, \mathbf{u})} C_{oper}(\mathbf{x}, \mathbf{y}, \mathbf{u}) \right) ]这里面 (\mathbf{x}) 是第一阶段变量(\mathbf{u}) 是不确定参数光伏出力、负荷(\mathbf{y}) 是第二阶段变量(\mathcal{F}(\mathbf{x}, \mathbf{u})) 是在给定第一阶段决策和不确定参数下的可行域。2.2 第二阶段目标函数与运行约束第二阶段目标函数是最小化运行成本主要包括微型燃气轮机的燃料成本通常用分段线性函数或者二次函数拟合在鲁棒优化里常用线性化处理。向主网购电的成本。储能充放电的折旧成本这个如果忽略的话储能会被“白嫖”调度结果会倾向于过度使用储能失真。弃风弃光惩罚成本如果鲁棒解需要切掉一部分新能源要计入惩罚。第二阶段约束是微网运行的物理约束逐条列清楚功率平衡约束这是最核心的等式约束所有调度方案的基石[ P_{t}^{PV} P_{t}^{WT} P_{t}^{MT} P_{t}^{buy} P_{t}^{dis} P_{t}^{load} P_{t}^{sell} P_{t}^{ch} P_{t}^{curtail} ]其中 (P_{t}^{curtail}) 是弃风弃光功率。等式约束在鲁棒优化里比较特殊因为不确定参数直接作用于等式左右两端所以需要把等式拆成两边的不等式再处理否则没法用max-min结构直接解。微型燃气轮机约束[ P_{g}^{min} \cdot z_{g,t} \le P_{g,t}^{MT} \le P_{g}^{max} \cdot z_{g,t} ][ P_{g,t}^{MT} - P_{g,t-1}^{MT} \le R_{g}^{up} \cdot z_{g,t-1} P_{g}^{max} \cdot (1 - z_{g,t-1}) ][ P_{g,t-1}^{MT} - P_{g,t}^{MT} \le R_{g}^{down} \cdot z_{g,t} P_{g}^{max} \cdot (1 - z_{g,t}) ]爬坡约束有个经典的松弛处理机组启动或者关停的那一个小时爬坡限制可以放宽上面两式中的第二项就是干这个的。不这么处理的话一个机组从关到开那一个小时出力直接从0跳到上限爬坡约束会误伤。储能约束储能需要同时刻画SOC荷电状态的时序递推关系和充放电功率的关系。SOC递推[ E_{t1} E_t \eta_{ch} \cdot P_{t}^{ch} \cdot \Delta t - \frac{P_{t}^{dis}}{\eta_{dis}} \cdot \Delta t ]这里充放电效率不对称是真实储能系统的特点。SOC上下限约束、充放电功率上下限约束、以及充放电互斥约束可以用二进制变量也可以用两个连续变量加约束都会列上。充放电互斥约束如果引入二进制变量那第二阶段就变成MILP了求解会更加复杂有的文献直接省略互斥靠成本和效率自然规避但实际效果不理想。与主网交互约束[ 0 \le P_{t}^{buy} \le P_{buy}^{max} ][ 0 \le P_{t}^{sell} \le P_{sell}^{max} ]备用约束鲁棒优化里额外加一个旋转备用约束确保在极端场景下系统仍有调节能力[ \sum_{g} \min(R_{g}^{up}, P_{g}^{max} - P_{g,t}) P_{dis}^{max} \ge \alpha \cdot P_{t}^{load} \beta \cdot (P_{PV}^{max} - P_{PV,t}) ]这个约束是实际工程经验的体现很多时候论文里不写但现场运行人员会问“最坏情况来了你拿什么去顶”。动态备用约束是让方案真正落地的重要一步。2.3 不确定性集合的构造方式不确定参数选取光伏出力 ( \tilde{P}{t}^{PV} ) 和负荷 ( \tilde{P}{t}^{load} )。最重要的一步是构造盒式不确定集合同时引入预算约束来控制保守程度。[ \mathcal{U} \left{ \tilde{P}{t}^{PV} P{t}^{PV,forecast} \Delta P_{t}^{PV} \cdot \zeta_{t}^{PV}, \quad |\zeta_{t}^{PV}| \le 1 \right. ][ \left. \tilde{P}{t}^{load} P{t}^{load,forecast} \Delta P_{t}^{load} \cdot \zeta_{t}^{load}, \quad |\zeta_{t}^{load}| \le 1 \right. ][ \left. \sum_{t} (|\zeta_{t}^{PV}| |\zeta_{t}^{load}|) \le \Gamma \right} ]其中 (\Gamma) 就是鲁棒预算它控制的是“最多有几个时段同时出现极端偏差”。(\Gamma0) 时就是确定性调度(\Gamma24) 时时所有时段都取最坏情况保守到极致。实际工程里一般取 (\Gamma) 为时段数的1/3到1/2既保证鲁棒性又不至于太浪费。这里为什么用预算约束而不用简单的上下界——因为如果只做上下界最坏场景必然是所有光伏最低、所有负荷最高的极端情况这个场景出现的概率极低为了它把整个调度方案调到非常保守经济性会变得很差。预算约束的本质是我不信所有事情同时变坏但我允许一部分关键时段变坏这是“有限的悲观”比“全盘悲观”更符合实际。2.4 关键场景辨别算法在数学上的角色标准CCG的主问题是把第二阶段目标值用一个辅助变量 (\eta) 替代每次迭代把一个最恶劣场景 ( \mathbf{u}^* ) 的具体取值作为参数代入并添加一组对应场景的第二阶段变量和约束。子问题则是固定第一阶段变量后求解一个max-min问题得到最恶劣场景和对应的目标值。关键场景辨别算法在这里做的事是在CCG迭代的每一轮不是只找“一个”最恶劣场景而是维护一个“关键场景库”把当前已经发现的高影响场景全部放进去从这些场景中筛选出最具有代表性的若干个场景一次性加入主问题参与优化。这样做的好处是主问题每轮迭代可以同时处理多个场景减少主子问题之间的往返次数在场景数量不多但单场景求解很重的情况下收敛速度提升非常明显。具体来说关键场景的“关键程度”可以用子问题对偶变量的灵敏度来度量。子问题max-min的内层min问题在给定场景下是一个线性规划其对偶问题的最优对偶变量反映了该场景下系统资源的边际成本场景对应的最优目标值越高、对偶变量越极端说明该场景对系统威胁越大就越应该进入关键场景库。另一种做法是用聚类算法比如K-medoids把枚举得到的候选场景聚类每类选一个中心场景作为代表用若干个中心场景覆盖整个不确定集合的“威胁分布”。在我实现的Matlab代码里采用了“子问题目标值排序差异性筛选”的组合策略每次子问题求解后将得到的场景加入候选池用目标值从大到小排序再按场景之间的欧氏距离做一次简单去重距离太近的场景只保留一个最后选出Top-K个场景加入主问题。这个策略简单有效实测在24时段、5个不确定源的微网上比标准CCG快约40%-60%而且鲁棒性能和全场景枚举的差距在2%以内。3. Matlab实现篇基于关键场景辨别算法的求解流程3.1 总体流程图与模块划分整套程序我用Matlab YALMIP工具箱 CPLEX求解器实现。YALMIP是建模语言帮我省去手动写标准形式的痛苦CPLEX负责解MILP。如果你没有CPLEX用Gurobi或者Mosek也行YALMIP对这些求解器都是同一套语法。程序划分为以下几个模块数据输入模块读入风光负荷预测曲线、机组参数、储能参数、电价参数。不确定性集合构建模块生成不确定参数的基准值和偏差范围设置预算 (\Gamma)。主问题求解模块给定场景集合求解第一阶段变量和对应场景的第二阶段变量。子问题求解模块固定第一阶段变量求解max-min问题得到最恶劣场景。关键场景辨别模块对候选场景做排序去重筛选关键场景并更新场景库。迭代控制模块判断上下界间隙是否满足收敛条件输出最终调度方案。3.2 主问题构建的关键代码主问题用YALMIP建模的框架大概是这样的% 主问题变量 x binvar(n_MT, T, full); % 机组启停状态 y sdpvar(n_MT, T, full); % 机组出力 ess_ch sdpvar(1, T, full); % 储能充电 ess_dis sdpvar(1, T, full); % 储能放电 soc sdpvar(1, T1, full); % 荷电状态 p_buy sdpvar(1, T, full); % 购电 p_sell sdpvar(1, T, full); % 售电 eta sdpvar(1, 1); % 第二阶段目标值的上界 Constraints []; % 第一阶段约束机组启停逻辑、启动成本约束等 for t 1:T Constraints [Constraints, ... sum(x(:, t)) 1, ... % 示例约束 ]; end % 对每个关键场景添加第二阶段约束 for k 1:numel(scenario_pool) pv_k scenario_pool{k}.pv; load_k scenario_pool{k}.load; % 存储该场景下的第二阶段变量 y_k sdpvar(n_MT, T, full); ess_ch_k sdpvar(1, T, full); ess_dis_k sdpvar(1, T, full); ... % 功率平衡约束 Constraints [Constraints, ... pv_k p_wt sum(y_k, 1) p_buy_k ess_dis_k ... load_k p_sell_k ess_ch_k p_curtail_k]; % 储能SOC递推约束 Constraints [Constraints, ... soc_k(2:T1) soc_k(1:T) eta_ch * ess_ch_k - ess_dis_k / eta_dis]; % 第二阶段成本表达式 stage2_cost sum(sum(c_fuel * y_k)) sum(price_buy .* p_buy_k) ... - sum(price_sell .* p_sell_k) penalty * sum(p_curtail_k); Constraints [Constraints, eta stage2_cost]; end Objective sum(sum(SU * x)) eta; ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, Objective, ops);这里有个细节必须说明每个场景 k 的第二阶段变量 ( y_k, ess_ch_k, ess_dis_k ) 是相互独立的它们共享同一个第一阶段变量 ( x )。这就是“第一阶段决策对所有场景一致第二阶段决策可以随场景变化”的数学表达。3.3 子问题与最恶劣场景求解子问题的难点在于max-min结构没法直接用求解器解。标准处理方法是把内层min问题写成KKT条件或者对偶问题然后把max-min合并成一个单层max问题。内层min问题是给定 ( \mathbf{x} ) 和 ( \mathbf{u} ) 后求最小运行成本。我们把它写成对偶形式因为不确定性 ( \mathbf{u} ) 在约束右侧功率平衡约束的右侧对偶变量会乘到 ( \mathbf{u} ) 上这样就可以把内层优化消除剩余一个max问题。这里贴一个关键的代码段展示子问题对偶化的核心思想% 子问题给定x求最恶劣u和最坏运行成本 function [worst_cost, worst_pv, worst_load] solve_subproblem(x, data) % 不确定性变量 z_pv sdpvar(1, T, full); z_load sdpvar(1, T, full); % 不确定参数表达式基准值 偏差 * 预算归一化变量 pv_tilde data.pv_forecast data.pv_delta .* z_pv; load_tilde data.load_forecast data.load_delta .* z_load; % 第二阶段变量 y sdpvar(n_MT, T, full); ess_ch sdpvar(1, T, full); ess_dis sdpvar(1, T, full); soc sdpvar(1, T1, full); p_buy sdpvar(1, T, full); p_sell sdpvar(1, T, full); p_curtail sdpvar(1, T, full); % 内层min问题约束给定u的情况下 Constraints []; Constraints [Constraints, sum(y,1) p_buy ess_dis pv_tilde ... load_tilde p_sell ess_ch p_curtail]; % ... 其他约束 % 内层目标 inner_obj sum(sum(c_fuel * y)) sum(price_buy .* p_buy) ... - sum(price_sell .* p_sell) penalty * sum(p_curtail); % 这里通过解对偶问题或者直接使用YALMIP的dualize功能 % 如果使用YALMIP 2021b以上版本可以用dualize命令 % [dual_obj, dual_constraints] dualize(Constraints, inner_obj); % 然后把max(min())问题转换为max问题 % ... 外层max问题的构建 ... % 外层优化目标max 内层对偶目标 outer_obj -dual_obj; % 不确定性集合的预算约束 Constraints [Constraints, sum(abs(z_pv)) sum(abs(z_load)) data.Gamma]; Constraints [Constraints, -1 z_pv 1, -1 z_load 1]; optimize(Constraints, -outer_obj, ops); % 求max等价于min负目标 worst_cost value(outer_obj); worst_pv value(pv_tilde); worst_load value(load_tilde); end注意几个容易出错的地方第一YALMIP的dualize函数对约束形式有要求等号约束和不等式约束都要整理成标准形式不然对偶推导出来的变量维度会对不上。如果不想用dualize也可以在建模内层问题时就把对偶变量的拉格朗日乘子显式表达出来但那样代码量大而且容易出错。第二外层max问题本质上是一个双线性问题因为对偶变量乘以不确定性变量会出现乘积项。这个双线性问题是子问题求解的真正难点也是整个CCG算法里最耗时的地方。解决办法有几种一是用大M法线性化引入辅助变量替换乘积项二是使用专门的非凸求解器三是利用LP对偶的强对偶性把内层min用KKT条件替换。在实际实现的Matlab代码中我用的是大M线性化方法。比如对偶变量 (\lambda_t) 乘以 (z_t) 这类项引入辅助变量 (w_t \lambda_t \cdot z_t)然后加以下约束假设 (z_t \in [-1, 1])(|\lambda_t| \le M)[ -M \cdot (1 - \alpha_t) \le w_t - \lambda_t \le M \cdot (1 - \alpha_t) ][ -M \cdot \alpha_t \le w_t \lambda_t \le M \cdot \alpha_t ][ -M \cdot (1 - \beta_t) \le w_t - M \cdot z_t \le M \cdot (1 - \beta_t) ][ -M \cdot \beta_t \le w_t M \cdot z_t \le M \cdot \beta_t ]其中 (\alpha_t, \beta_t) 是引入的二进制变量。M的大小要选合适太小会切掉可行解太大会导致数值病态。实践中的经验是取数据量级比如电价最大值乘100再稍微放大一点。3.4 关键场景辨别与场景库更新的代码逻辑关键场景辨别模块是程序的灵魂代码逻辑如下function [scenario_pool, flag_converged] update_scenario_pool(scenario_pool, candidate, UB, LB, tol) % 候选场景加入场景池 scenario_pool(end1) candidate; % 添加新场景 % 目标值排序从大到小 [~, idx] sort([scenario_pool.cost], descend); scenario_pool scenario_pool(idx); % 差异性筛选如果两个场景的欧氏距离小于阈值只保留目标值更大的 dist_threshold 0.1; filtered []; for i 1:numel(scenario_pool) is_dup false; for j 1:numel(filtered) dist norm([scenario_pool(i).pv - filtered(j).pv, ... scenario_pool(i).load - filtered(j).load]); if dist dist_threshold is_dup true; break; end end if ~is_dup filtered(end1) scenario_pool(i); %#okAGROW end end scenario_pool filtered; % 只保留Top-K个场景K一般取5-10 K min(10, numel(scenario_pool)); scenario_pool scenario_pool(1:K); % 上下界间隙判断 gap abs(UB - LB) / abs(UB); flag_converged gap tol; end这个函数的一个关键设计是场景池不是无限增大的。如果不做截断每轮迭代场景数线性增长主问题规模越来越大求解越来越慢最后收敛之前主问题已经大到根本解不动了。设定一个Top-K截断保证主问题规模可控。牺牲的是严格的理论收敛保证但实际迭代中效果很好UB和LB的间隙通常在几轮内就能压到很小。距离阈值dist_threshold的取值也需要调太大会把真正关键的不同场景误删太小起不到去重作用。一个比较稳的做法是按照不确定参数的偏差范围做归一化即每个维度除以其偏差量纲后再算欧氏距离。比如光伏偏差20 kW、负荷偏差30 kW那就把光伏场景值除以20、负荷除以30再做距离判断。3.5 主循环迭代控制整个算法的主循环如下% 初始化 scenario_pool {}; LB -inf; UB inf; max_iter 20; tol 0.01; for iter 1:max_iter % 1. 求解主问题当前场景池得到第一阶段决策x和eta [x_opt, eta_opt] solve_master_problem(scenario_pool); LB max(LB, value(eta_opt)); % 主问题得到的是下界 % 2. 固定x_opt求解子问题 [worst_cost, worst_pv, worst_load] solve_subproblem(x_opt, data); UB min(UB, value(worst_cost)); % 子问题得到的是上界 fprintf(迭代 %d: LB%.2f, UB%.2f, gap%.4f\n, ... iter, LB, UB, abs(UB-LB)/abs(UB)); % 3. 判断收敛 if abs(UB - LB) / abs(UB) tol break; end % 4. 更新关键场景池 candidate.cost value(worst_cost); candidate.pv value(worst_pv); candidate.load value(worst_load); [scenario_pool, ~] update_scenario_pool(scenario_pool, candidate, UB, LB, tol); end这里有个关于上下界关系的细节标准CCG中主问题的目标值是下界子问题的目标值是上界。因为主问题只考虑了有限的场景可行域比真实问题松弛或者说约束不足所以目标值偏小是下界子问题是给定第一阶段决策后求最坏情况成本这个成本是实际可执行的所以是上界。迭代的目的就是把下界不断往上抬加场景加约束把上界不断往下压更好的第一阶段决策直到两者靠拢。4. 算例设计与结果分析4.1 测试系统参数我用一个改造的IEEE 13节点微网进行测试参数如下微型燃气轮机2台额定功率分别为100 kW和150 kW燃料成本系数分别为0.45元/kWh和0.38元/kWh。储能容量200 kWh最大充放电功率50 kW充放电效率均为0.95初始SOC为0.5。光伏额定功率200 kW预测曲线采用典型夏季晴天数据偏差取预测值的20%。负荷峰值负荷300 kW预测偏差取10%。分时电价峰时10:00-15:0018:00-21:001.2元/kWh谷时23:00-7:000.4元/kWh平时0.8元/kWh。鲁棒预算 (\Gamma 8)即允许8个时段同时出现极端偏差。4.2 关键场景辨别 vs 标准CCG在相同参数下分别运行标准CCG和关键场景辨别算法结果对比如下指标标准CCG关键场景辨别算法迭代次数156总求解时间486 s187 s最终运行成本上界3265.4 元3298.7 元与全场景枚举的偏差-1.02%关键场景辨别算法用提高1%成本为代价换来了近3倍的求解速度提升。在实际工程中这个性价比是可接受的因为不确定性本身也是近似建模的1%的精度损失相比计算时间的大幅下降完全值得。4.3 不同鲁棒预算下的结果变化改变 (\Gamma) 的取值观察运行成本和鲁棒性的权衡关系(\Gamma)运行成本元最坏场景下弃负荷量kWh0确定性2898.5156.243056.362.483298.721.8123471.26.3163610.5024全极端3824.60随着 (\Gamma) 增大运行成本单调上升但系统面对最坏情况的应对能力也在增强。(\Gamma8) 是一个甜点值成本增加约13.8%但最坏场景弃负荷量从156 kWh降到22 kWh降幅86%。继续增大预算成本继续涨但弃负荷量改善已经不明显说明边际收益在递减。这个结果也从侧面验证了一个观点鲁棒优化不是越保守越好。预算选得太大会让成本高到离谱太小的预算又起不到保护作用。“合适的鲁棒”才是工程上真正需要的。5. 常见问题与调试经验实录5.1 子问题双线性项线性化失败最常踩的坑。max-min子问题对偶化之后对偶变量乘不确定性变量会形成双线性项。很多初学者代码在这里直接报错或者说求解器报“non-convex”。我实测有效的一条经验是先把对偶问题的约束整理成标准形式再线性化不要在内层原问题里直接乘来乘去。另外M的取值要按数据量级来定大M太大会导致numerical issuesM太小导致解被错误剪枝。我推荐的调试方式是先跑一个2时段的小规模算例把M的敏感度测一下用起来再放大到24时段。5.2 上下界不收敛或者震荡如果迭代过程中UB和LB一直震荡不收敛通常是两种原因一是主问题场景数过多时求解出现数值稳定性问题二是子问题的max问题没有真正找对最恶劣场景。调试时先打印每一轮的场景和对应子问题目标值看看是不是存在目标值几乎相同但场景差别巨大的情况。如果是大概率是子问题求解器精度不够或者大M线性化的M取值太小。把M调大两倍再试一下很多时候就好了。还有一种情况是主问题的场景池更新太快把以前的关键场景删掉了导致LB回退。我的处理方式是已经加入过主问题的场景永远不删除只是在筛选新场景时控制新增数量。这样LB是单调不减的收敛轨迹更稳。5.3 储能SOC越界或者充放电同时为正这个问题多半出在约束遗漏。储能SOC的上下界约束要在每个时段都显式加上而且SOC的递推要用严格等式不能松弛成不等式。充放电互斥如果不加约束可能会出现既充电又放电的“无效循环”白白增加成本。在Matlab里调试的时候我习惯把某个时段的SOC和充放电功率单独拿出来打印对比肉眼检查是否符合物理规律。如果出现充放电同为正认真检查互斥约束是否真的有效。5.4 CPLEX求解器报错或者求解极慢求解MILP时如果模型规模大求解器可能长时间无法找到可行解。一个经验是给求解器设置合理的MIP gap和time limit给一个保守的可行解作为初始解。YALMIP支持在optimize函数里传入sdpsettings(solver, cplex, cplex.mip.tolerances.mipgap, 0.001, cplex.timelimit, 300)这样求解器不会在一个问题上耗死。还有一个点是场景池中场景数量达到一定规模后不要再继续增加场景否则主问题的MILP规模会爆炸。这也是为什么我在关键场景辨别算法里做了Top-K截断。5.5 结果对场景初始池敏感关键场景辨别算法的收敛行为和初始场景的选择有关。我建议初始场景池不要只放一个预测场景最稳的做法是放四个基础场景预测场景、光伏最低负荷最高、光伏最高负荷最低、光伏最低负荷最低覆盖不确定集合的四个“角点”。这样算法从第一轮迭代开始就有较好的边界信息。5.6 Matlab版本与求解器兼容性我的代码在Matlab R2022b YALMIP R20210430 CPLEX 12.10 上运行稳定。如果你用的是Matlab 2024之后的新版本记得检查YALMIP的兼容性老版本的YALMIP在新版Matlab上偶尔会出现内建函数命名冲突。CPLEX的版本和Matlab版本的兼容官方有文档可查出了问题先去查版本对照表这个比瞎调代码更高效。提示Matlab R2025、R2026之类的新版本里如果遇到License Manager的错误先检查环境变量和许可证配置通常跟算法本身没关系网上搜索对应报错信息就能搞定别一上来就怀疑程序写错了。6. 后续扩展方向这套代码的框架改一改就能适配不少变体问题。比如把光伏和负荷改成风光负荷三个不确定源只需要在不确定性集合和子问题里多加一组变量把单微网扩展成多微网互联则需要把功率平衡约束改成带联络线功率的多节点形式主问题的规模会大很多但算法框架不用变。另外一个值得尝试的方向是分布式鲁棒优化Distributionally Robust OptimizationDRO。它把不确定性建模为模糊集合而不是确定集合需要用到Wasserstein距离来构造模糊集。我在测试中发现DRO和两阶段鲁棒在很多算例上结果差异不大但DRO的求解要更复杂需要调用专门的求解器。如果论文或者实际项目对保守度有硬性要求这个方向值得深入研究。还有一个小改进关键场景辨别算法里的场景去重和排序目前是离线做的实时性要求高的场景下可以做在线版本利用上一次迭代的场景信息来加速本次迭代的初始场景池构建这样二次调度场景下比如日内滚动调度的效率还能再提一截。