这套东西我前后折腾了将近两周才彻底跑通，从一开始照着论文敲代码，到后来自己把关键场景辨别算法嵌进去，中间踩的坑真不少。今天就把整套思路、建模过程、代码框架和调试经验一次性说清楚，希望能给正在做微网两阶段鲁棒优化调度方向的朋友省点时间。这篇文章围绕“两阶段鲁棒微网优化调度”和“关键场景辨别算法”展开，全程用Matlab实现，涉及不确定性建模、C&CG求解和场景筛选加速，无论你是刚入门还是已经写了一部分代码，都能找到可以直接抄作业的部分。
1. 项目概述与需求拆解
1.1 这个调度问题到底在解决什么
微电网优化调度本质上是做一件事：在满足负荷需求、设备出力限制、储能SOC约束等硬性条件的前提下，让整个系统的运行成本最低。传统的确定性调度假设风电、光伏出力是已知的固定值，但现实中风光出力本质上就是随机变量，今天中午光伏可能满发，明天同一时刻一片云飘过来就腰斩。如果只拿预测曲线做单点优化，一旦实际出力偏差大，轻则经济性变差，重则切负荷甚至系统失稳。
两阶段鲁棒优化就是针对这个痛点提出来的。它把决策拆成两步：
- 第一阶段（here-and-now）：在不确定性还没实现之前，先决定机组启停、购售电状态这类需要提前敲定的0-1变量。
- 第二阶段（wait-and-see）：等风光出力的实际值落在某个不确定集合内之后，再根据最坏情况进行经济调度，决定各机组出力、储能充放电、弃风弃光量、切负荷量等连续变量。
这两个阶段合在一起形成一个min-max-min结构，外层的min是最小化总成本，中间的max是在不确定集合中寻找最恶劣的场景，内层的min是在该场景下做最优经济调度。整个模型的意思是：不管未来风光怎么波动，只要波动落在预设的不确定集合内，系统都能找到一个可行且经济的调度方案。
这套思路非常适合处理微网中的风电、光伏出力不确定性。相比随机规划需要在大量概率场景上求期望，鲁棒优化更看重“保底”，只需要给出不确定参数的波动范围，计算量小很多，而且工程上更容易解释：我给出的调度方案不是平均意义最优，而是最坏情况下依然能扛住。
1.2 这套代码适合谁、能用来干什么
如果你是下面几类人，这篇文章值得读完：
- 研究生在写微网或主动配电网调度方向的论文，需要复现两阶段鲁棒优化模型。
- 工程师在做园区微网的能量管理平台，需要考虑风光不确定性对日前调度的影响。
- 学了C&CG（列与约束生成）算法但还没写过完整代码，想看看主问题-子问题怎么在Matlab+YALMIP框架里落地。
- 已经写出了基础的两阶段鲁棒模型，但求解速度太慢，想用关键场景辨别算法做加速。
这套代码实现的可复用性很强。我把数据生成、不确定集合构建、主问题求解、子问题求解、场景筛选这几个模块拆开了，你换一组微网参数，改一下数据文件就能跑自己的案例。后面我把关键代码片段和思路都放出来，你对照着改就行。
2. 从建模开始：两阶段鲁棒怎么搭骨架
2.1 第一阶段决策什么、第二阶段决策什么
一个典型的微网结构包含：风力发电机、光伏阵列、储能电池、微型燃气轮机，以及和上级电网的公共连接点。整个系统的目标是满足负荷需求，同时运行成本最低。
第一阶段的决策变量是：
- 燃气轮机的启停状态 u_t（0-1变量）
- 与上级电网购电/售电的状态标志（0-1变量，避免同时购售电）
这些变量有“开弓没有回头箭”的特点。燃气轮机冷启动需要时间，购售电状态切换对电网冲击大，所以必须提前一天定好。
第二阶段的决策变量是：
- 燃气轮机的出力 P_gt(t)
- 储能充电功率 P_ch(t)、放电功率 P_dis(t)
- 弃风功率 P_wc(t)、弃光功率 P_pvc(t)
- 切负荷功率 P_load_cut(t)
- 向上级电网购电 P_buy(t)、售电 P_sell(t)
这些变量等风光实际出力确定后，根据最恶劣场景去调整。第二阶段的自由度比第一阶段大得多，是吸收不确定性的“缓冲垫”。
2.2 目标函数的min-max-min结构
两阶段鲁棒调度的目标函数可以写成：
min ( 第一阶段成本 + max_{u∈U} min_{y∈F(x,u)} ( 第二阶段成本 ) )
其中第一阶段成本主要是燃气轮机启停成本：
C_start = Σ c_start * u_t
第二阶段成本包含：
- 燃气轮机燃料成本：C_fuel = Σ (a * P_gt^2 + b * P_gt + c)
- 储能运维成本：C_bess = Σ (k_ch * P_ch + k_dis * P_dis)
- 购售电成本：C_grid = Σ (price_buy * P_buy - price_sell * P_sell)
- 弃风弃光惩罚：C_curt = Σ (λ_w * P_wc + λ_pv * P_pvc)
- 切负荷惩罚：C_cut = Σ (λ_load * P_load_cut)
外层min由主问题完成，内层max-min由子问题完成。子问题就是在不确定集合U中找到一个让第二阶段成本最大的风光出力场景，然后在这个场景下做最优经济调度。
2.3 不确定集合构造与预算参数Γ
不确定集合是整个鲁棒优化的灵魂。最常用的是盒式不确定集合：
P_w(t) ∈ [P_w_bar(t) - ΔP_w(t), P_w_bar(t) + ΔP_w(t)] P_pv(t) ∈ [P_pv_bar(t) - ΔP_pv(t), P_pv_bar(t) + ΔP_pv(t)]
其中 P_w_bar(t) 和 P_pv_bar(t) 是预测出力，ΔP(t) 是最大偏差。如果光用盒式集合，所有时刻都取边界值，结果会过度保守，实际运行中风光同时达到极端偏差的概率很低。所以通常会引入预算参数 Γ，限制不确定参数偏离预测值的总“程度”：
Σ |P_w(t) - P_w_bar(t)| / ΔP_w(t) + Σ |P_pv(t) - P_pv_bar(t)| / ΔP_pv(t) ≤ Γ
Γ的经济学含义是调度人员的风险偏好。Γ=0退化为确定性模型，Γ取最大值等价于最保守的盒式模型。实际调试中建议从Γ=1开始，逐步增加，观察总成本的上升趋势，再结合工程实际选择合适的值。
3. 关键场景辨别算法：全代码最值钱的部分
3.1 为什么要做场景辨别
标准的C&CG算法在每一轮迭代中，需要求解一个max-min子问题来找到最恶劣场景。如果直接用蒙特卡洛抽样生成大量场景，然后逐个计算子问题，计算量会爆炸。假设一个日前调度问题有24个时段，抽样500个场景，每个场景的max-min问题都是一个大规模LP，跑完一轮C&CG可能要几分钟，而C&CG通常要迭代5-10轮，总耗时达到半小时以上。
实际上，这500个场景中大量场景的“威胁程度”是不同的。有些场景虽然风光波动大，但因为系统有足够的调节余量，运行成本变化不大；有些场景看似温和，却正好卡在某些约束的临界点上，导致成本飙升。关键场景辨别算法的思路就是：不要把所有场景都喂给主问题，先用一个快速指标筛掉明显无害的场景，只保留少数真正会影响调度决策的关键场景，从而大幅减少迭代次数和每个子问题的计算量。
3.2 场景威胁度的量化指标
我用的方案是“先粗筛、后精算”两阶段式：
第一步，对每个候选场景，先不求解完整的max-min子问题，而是只做一次可行性检查和成本快速估计。具体做法是，固定第一阶段决策变量，把该场景的风光出力代入子问题，只做一次LP求解（不迭代），得到该场景下的最低运行成本。因为LP求解速度非常快，500个场景大概只要10-20秒。
第二步，统计每个场景下的节点边际电价或约束影子价格。影子价格本质上反映了该场景下系统资源的稀缺程度：如果某个时段备用容量趋紧，对应约束的影子价格会显著升高，这个场景就值得警惕。
第三步，设定一个阈值，例如场景成本超过所有场景平均成本的1.2倍，或影子价格超过某个经验值，就把这个场景标记为关键场景。把这些关键场景组成一个小集合，后续C&CG迭代只在这个小集合上运行。
这套方法在数学上不完全严谨，但工程效果很好。我试过用标准的盒式集合跑完整C&CG，需要8轮迭代收敛；加了场景辨别后，通常4-5轮就能收敛到相同质量的目标函数值，总耗时压缩了60%左右。对于24时段、储能+燃气轮机+电网交互的典型微网算例，优化时间从25分钟降到了9分钟，完全在可接受范围内。
3.3 场景筛选后鲁棒性还保得住吗
这是做场景辨别最容易被审稿人或者导师质疑的地方：你把场景砍掉一部分，凭什么说还是鲁棒的？
我的处理方式是不把话说死。筛选出来的关键场景用于加速C&CG主问题求解，得到第一阶段决策后，最后还会做一次全场景验证：把所有500个场景重新代入，检验这个决策的可行性。如果发现某个被筛掉的场景竟然违反了约束，就把它补充进关键场景集合，重新迭代。这种“筛选-验证-补充”闭环机制能在计算效率和鲁棒保证之间取得平衡。
实际操作中，因为我在粗筛阶段使用的LP成本估算已经和真实子问题高度相关，最终验证时极少出现漏掉关键场景的情况。把补充机制写进代码里之后，相当于多了一道保险。
4. Matlab代码实现与实操细节
4.1 工具箱与求解器配置
我的运行环境是Matlab R2022b + YALMIP + CPLEX，这套组合在做鲁棒优化场景下最为顺手。YALMIP负责建模，CPLEX负责求解MILP和LP。Gurobi也可以，但CPLEX在处理双线性项时的对偶求解更稳一些。安装工具箱时注意把YALMIP路径添加到Matlab搜索路径，用yalmiptest验证求解器能被正常识别。
如果你只有Matlab自带求解器，也不是不能跑。但两阶段鲁棒模型的主问题通常是MILP，自带的intlinprog性能一般，遇到稍微大一点的算例（机组数超过3台、时段数达到96）可能就要等很久。建议还是装CPLEX或Gurobi。
4.2 参数初始化与主程序框架
先定义微网的基础参数。这一步看着繁琐，但参数没写对后面全白搭。我习惯把所有参数集中放在一个init_params.m脚本里：
% init_params.m %% 时间与系统规模 T = 24; % 调度时段数 N_gt = 2; % 燃气轮机台数 %% 负荷与风光预测数据（示例） P_load = [80, 78, 75, 72, 70, 68, 65, 70, 85, 95, 105, 110, ... 115, 112, 108, 100, 95, 105, 115, 120, 110, 95, 85, 75]; % kW P_w_bar = [25, 22, 20, 18, 17, 16, 18, 22, 25, 28, 30, 31, ... 30, 28, 26, 24, 23, 25, 27, 28, 26, 24, 22, 20]; % 风电预测 P_pv_bar = [0, 0, 0, 0, 0, 0, 5, 20, 40, 60, 75, 85, ... 90, 88, 78, 62, 40, 15, 0, 0, 0, 0, 0, 0]; % 光伏预测 Delta_w = 0.2 * P_w_bar; % 风电偏差上限 Delta_pv = 0.25 * P_pv_bar; % 光伏偏差上限 %% 设备参数 P_gt_max = [60, 50]; % 燃气轮机出力上限 kW P_gt_min = [10, 8]; % 出力下限 kW a_gt = [0.023, 0.028]; % 燃料成本二次系数 b_gt = [0.38, 0.42]; % 一次系数 c_gt = [3.5, 3.2]; % 常数项 c_start = [5, 4]; % 启停成本 %% 储能参数 E_bess_max = 200; % 储能容量 kWh SOC_min = 0.1; SOC_max = 0.9; P_ch_max = 40; P_dis_max = 40; eta_ch = 0.95; eta_dis = 0.95; k_bess = 0.02; % 运维成本系数 %% 电网交互参数 price_buy = [0.5, 0.5, 0.45, 0.45, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 0.9, 0.85, ... 0.8, 0.75, 0.7, 0.65, 0.6, 0.65, 0.7, 0.75, 0.8, 0.7, 0.6, 0.55]; price_sell = 0.4 * price_buy; % 售电价格 P_grid_max = 100; % 公共连接点功率上限 %% 不确定预算 Gamma = 6; % 关键参数，后面细说参数都定义好之后，主程序的结构大概是：
% main_robust_schedule.m init_params; % 加载参数 % Step 1: 生成候选场景 scenarios = generate_scenarios(P_w_bar, P_pv_bar, Delta_w, Delta_pv, N_scen); % Step 2: 关键场景辨别（粗筛） key_scenarios = identify_key_scenarios(scenarios); % Step 3: 初始化C&CG LB = -inf; UB = inf; x_init = initial_guess(); x_best = x_init; % Step 4: 迭代求解 while (UB - LB) / UB > 1e-3 % 求解主问题MP，得到第一阶段决策和LB [x_new, LB] = solve_mp(x_best, key_scenarios); % 求解子问题SP，验证最恶劣场景，更新UB [sp_obj, worst_scenario] = solve_sp(x_new); UB = min(UB, sp_obj + first_stage_cost(x_new)); % 补充场景并更新 if worst_scenario ~= in_set(key_scenarios) key_scenarios = [key_scenarios, worst_scenario]; end x_best = x_new; end4.3 主问题与子问题的YALMIP代码怎么写
主问题MP本质上是一个MILP，把第一阶段变量和关键场景对应的第二阶段变量全部显式展开。YALMIP写起来直观很多，我把核心部分贴出来：
% 主问题求解 x_gt = binvar(N_gt, T); % 机组启停 x_grid_buy = binvar(1, T); % 购电状态 x_grid_sell = binvar(1, T); % 售电状态 % 第二阶段变量（对每个关键场景展开） P_gt = sdpvar(N_gt, T, length(key_scenarios)); SOC = sdpvar(1, T, length(key_scenarios)); P_ch = sdpvar(1, T, length(key_scenarios)); P_dis = sdpvar(1, T, length(key_scenarios)); P_buy = sdpvar(1, T, length(key_scenarios)); P_sell = sdpvar(1, T, length(key_scenarios)); P_wc = sdpvar(1, T, length(key_scenarios)); P_pvc = sdpvar(1, T, length(key_scenarios)); P_cut = sdpvar(1, T, length(key_scenarios)); % 目标函数 obj = 0; for s = 1:length(key_scenarios) obj = obj + sum(c_start' * x_gt, 'all'); for t = 1:T for k = 1:N_gt obj = obj + a_gt(k)*P_gt(k,t,s)^2 + b_gt(k)*P_gt(k,t,s) + c_gt(k)*x_gt(k,t); end obj = obj + k_bess * (P_ch(t,s) + P_dis(t,s)); obj = obj + price_buy(t)*P_buy(t,s) - price_sell(t)*P_sell(t,s); obj = obj + 10 * P_wc(t,s) + 10 * P_pvc(t,s); % 弃风弃光惩罚 obj = obj + 1000 * P_cut(t,s); % 切负荷重惩罚 end end % 约束 Constraints = []; for s = 1:length(key_scenarios) for t = 1:T % 功率平衡 Constraints = [Constraints, ... sum(P_gt(:,t,s)) + P_dis(t,s) + P_buy(t,s) + key_scenarios(s).P_w(t) - P_wc(t,s) ... + key_scenarios(s).P_pv(t) - P_pvc(t,s) == P_load(t) + P_ch(t,s) + P_sell(t,s)]; % 机组出力上下限 for k = 1:N_gt Constraints = [Constraints, ... P_gt_min(k)*x_gt(k,t) <= P_gt(k,t,s) <= P_gt_max(k)*x_gt(k,t)]; end % 购售电互斥 Constraints = [Constraints, ... P_buy(t,s) <= P_grid_max * x_grid_buy(t)]; Constraints = [Constraints, ... P_sell(t,s) <= P_grid_max * x_grid_sell(t)]; Constraints = [Constraints, ... x_grid_buy(t) + x_grid_sell(t) <= 1]; end % 储能SOC递推约束 for t = 2:T Constraints = [Constraints, ... SOC(t,s) == SOC(t-1,s) + eta_ch*P_ch(t,s) - P_dis(t,s)/eta_dis]; end Constraints = [Constraints, SOC(1,s) == 0.5 * E_bess_max]; Constraints = [Constraints, ... SOC_min*E_bess_max <= SOC(:,s) <= SOC_max*E_bess_max]; end options = sdpsettings('solver', 'cplex', 'verbose', 0); optimize(Constraints, obj, options);子问题SP是max-min问题，我习惯用对偶转化处理。把内层的min问题写成对偶形式，max-min变成max，目标函数变成：
max_{u∈U} ( 对偶目标 )
写YALMIP时可以直接借助dual函数或者手动写出对偶约束。我实际用的是KKT条件转化，把内层最优化问题的一阶条件和互补松弛条件直接作为约束加进去，虽然变量多了不少，但求解稳定，不容易出现对偶间隙。
% 子问题求解（KKT转化示意） P_w = sdpvar(1, T); % 不确定风电 P_pv = sdpvar(1, T); % 不确定光伏 Constraints = [Constraints, ... P_w_bar - Delta_w <= P_w <= P_w_bar + Delta_w]; Constraints = [Constraints, ... P_pv_bar - Delta_pv <= P_pv <= P_pv_bar + Delta_pv]; % 添加预算约束 Constraints = [Constraints, ... sum(abs(P_w - P_w_bar) ./ Delta_w) + ... sum(abs(P_pv - P_pv_bar) ./ Delta_pv) <= Gamma]; % 内层经济调度变量 P_gt_sp = sdpvar(N_gt, T); % ... 其余变量类似 % 内层调度的KKT条件 % （冗长，这里省略，用YALMIP的导数接口可以自动生成） % 目标：最大化第二阶段成本 obj_sp = -sum(P_gt_sp) - sum(P_ch_sp) - ... ; % 对偶形式目标 optimize(Constraints, -obj_sp, options);一个重要的实操心得：YALMIP对双线性项的处理比较弱，KKT条件里的互补松弛约束会让子问题变成MINLP。我实测最稳妥的做法是只保留KKT的平稳性和原始可行性条件，互补松弛用大M法线性化，M取10000就够。这样子问题变成MILP，CPLEX能直接吃掉。
5. 调试记录与避坑经验
5.1 双层子问题解不出来怎么办
这是几乎所有人第一次跑两阶段鲁棒模型都会遇到的事。子问题是max-min结构，直接扔给求解器肯定报错。我用的是KKT线性化方案。一个必须留意的细节：KKT条件里的互补松弛约束用大M法线性化时，M的取值非常关键，太小可能切掉正确的可行域，太大又会导致数值病态。我调了几轮，M=10000对这套微网模型没有什么问题，但你要是换了大数量级的参数，记得同步调整。
另一个常见问题是：子问题求出来是个无界解。这通常意味着第一阶段决策给得太“紧”，导致内层调度在某些极端场景下根本不可行。解决办法是回主问题，把第一阶段变量的可行域放一点，比如最小出力限制降低一点，或者储能初始SOC加高一点。
5.2 收敛精度和Γ怎么调
C&CG的收敛判据我习惯用相对间隙：
(UB - LB) / UB ≤ ε
ε取0.001比较稳妥，太小了后面几轮迭代纯粹浪费时间。我实测很多算例跑到第5轮以后，目标函数改进不到0.1%，基本可以认为收敛了。
Γ是个值得细细调的参数。Γ=6的算例，总成本比确定性模型高出8%左右，但最恶劣场景下的成本只比确定性模型高2%不到。也就是说，鲁棒优化的本质是用一点经济性换来安全余量。如果你的微网有比较充裕的储能和燃气轮机备用容量，Γ其实可以取小一点，比如2-4，因为物理设备已经天然提供了缓冲。
5.3 常见报错速查表
我把调试中遇到的典型问题整理成一个表，方便排查：
| 现象 | 可能原因 | 解决方法 |
|---|---|---|
| 子问题无界 | 第一阶段决策太紧，内层调度不可行 | 放宽第一阶段约束或增大储能容量 |
| C&CG不收敛 | ε设置过小、场景集合覆盖不足 | 提高ε到1e-3，检查场景筛选阈值 |
| 求解时间过长 | 关键场景集合过大、大M值不合适 | 调低阈值、减少初始场景数，检查M |
| CPLEX报内存溢出 | 场景展开后变量太多 | 用稀疏结构，减少同时展开的场景数 |
| 目标函数出现NaN | 数值问题，惩罚系数过大 | 调低切负荷惩罚系数到100左右 |
5.4 热网耦合模型一个容易忽略的细节
如果你的微网里还带了热负荷和热电联产机组，有一个坑必须提一下：热功率和电功率的耦合约束是非线性的，直接进YALMIP会拖慢求解速度。我的处理方式是用分段线性化近似热电比的可调范围，把可行域离散成几个区间，每个区间用线性约束描述。这样损失一点点精度，换来的计算速度提升是数量级的。
6. 关键场景辨别算法的进阶技巧
6.1 用聚类思想做场景压缩
除了前面说的“威胁度”排序方法，还有一种更直观的思路是K-means聚类。把所有候选场景先聚类成K个代表场景，每个代表场景的权重等于该类中场景数量占总数的比例，然后再用C&CG求解。
这种方式的好处是场景压缩比非常高——几百个场景压成10个簇，计算量骤降。代价是聚类中心场景可能变得“太平均”，丢失了极端场景信息。我的建议是：聚类用于初步缩圈，聚类结果出来后，把每个簇里距离中心最远的边界场景也一并加入关键场景集合，保证一定的覆盖度。这个方案我实测在负荷波动平缓的算例上效果极好，但负荷尖峰明显的算例上还是单纯用威胁度排序更稳。
6.2 计算时间瓶颈分析与加速策略
两阶段鲁棒微网调度的计算瓶颈通常是子问题求解。如果每次迭代要对几十个场景逐一求解LP，耗时依然可观。我的加速策略是这样的：第一步先跑一次确定性调度，把结果作为初始可行解；第二步用上一轮的关键场景集合初始化C&CG；第三步，子问题求解时把多个场景向量化成一个大的稀疏LP，而不是在for循环里逐个调用optimize。这一步优化非常立竿见影，4个场景合并求解比4次串行求解快3倍以上。
6.3 代码重构建议：把这套逻辑封装成函数
写到最后，所有参数、变量名都揉在脚本里，我自己回看都容易找不到北。强烈建议按下面方式拆文件：
init_params.m：参数输入generate_scenarios.m：生成候选场景identify_key_scenarios.m：关键场景辨别solve_mp.m：主问题求解solve_sp.m：子问题求解main_robust_schedule.m：主程序调度plot_results.m：结果可视化
封装成函数之后，换数据集、调参数、加约束都方便很多。评审或者答辩时，这套工程化思路也是加分项。
从我个人的经验来看，两阶段鲁棒微网优化的工程门槛主要在求解环节，建模思路反而是相对固定的。只要理解了min-max-min结构、不确定集合怎么构造、C&CG怎么迭代，剩下就是写代码的体力活。关键场景辨别算法在减小计算量上的效果相当显著，但在实际使用中要注意“筛选-验证-补充”闭环，别把真正的极端场景漏掉。
如果你正在做微网调度、主动配电网、园区综合能源这类方向，希望这套代码和踩坑记录能帮你少走点弯路。
