两阶段鲁棒微网优化调度：关键场景辨别与C&CG算法实战

做微网优化调度这几年，最常被同行问的一个问题是：风光出力曲线预测不准，做出来的调度方案到底能不能扛住极端天气？确定性优化在预测值上做得再漂亮，一旦实际风光出力偏离个两三成，储能和燃气轮机的调节压力立刻拉满，严重的时候甚至得切负荷。这也是我后来把研究方向转向两阶段鲁棒微网优化调度的直接原因——与其赌预测精度，不如直接承认不确定性存在，然后用鲁棒优化的框架把最坏情况兜住。

这篇文章复盘的就是一个完整可跑的方案：基于关键场景辨别算法的两阶段鲁棒微网优化调度，代码用Matlab实现。内容覆盖了从不确定性建模、关键场景筛选、C&CG求解算法到代码模块拆解的完整链路。不管你是刚开始接触鲁棒优化的研究生，还是已经在做微网能量管理、想把手里的确定性调度模型升级成鲁棒版本的工程师，这篇都值得花二十分钟看完。

1. 两阶段鲁棒微网调度：这个课题到底在解决什么问题

1.1 为什么微网调度必须做成“两阶段”

先说清楚两阶段的意思。微网调度本质上是一个“决策-等待-再决策”的过程：在日前阶段，你要决定燃气轮机的开停机、是否需要预留备用容量、储能是否进入待命状态，这些是离散决策，定了之后短时间内改不了；到了日内运行阶段，风光出力、负荷逐渐变成已知量，你再根据实际情况调整燃气轮机出力、储能充放电功率、与上级电网的购售电功率。前者叫第一阶段决策（here-and-now），后者叫第二阶段决策（wait-and-see）。

如果只做一个单阶段优化，把预测值当成真实值，那相当于让调度员在不知道未来天气的情况下把所有决策一次性拍死。风光出力超预测时储能可能没空间充电，出力不及预期时燃气轮机又可能爬坡跟不上。两阶段结构的核心价值在于：第一阶段只做那些必须提前定的决策，把大量连续调节的灵活性留到第二阶段，不确定参数观测到之后再做适应性调整。这个结构天然贴合实际运行逻辑，也是微网鲁棒调度普遍采用两阶段框架的根本原因。

1.2 确定性调度的短板与不确定性问题的本质

确定性调度在工程上实现简单，Yalmip建模、Cplex一跑，几分钟就出结果。但它的前提是预测误差很小，调度方案对预测偏差的容忍度却很低。实操中经常出现这样的情况：某天光伏预测出力为100kW，确定性优化据此安排燃气轮机少发、储能少充，结果午后云层突然加厚，光伏实际出力只有60kW，储能因为前一天已经按计划充满而无处充电，燃气轮机爬坡上限又不够，最终只能向主网购电或者切负荷。

随机优化是另一条路，但它要求你准确知道各随机变量的概率分布，而微网层面的风光负荷历史数据往往不够充分，分布估计不准时优化结果反而误导运行决策。鲁棒优化则换了一种思路：我不需要精确分布，只需要知道每个不确定参数的可能波动范围，然后保证在这个范围内的任何实现下方案都可行。代价是结果偏保守，但这个保守是可以接受且可以调节的——通过预算Γ控制偏离程度，你可以在鲁棒性和经济性之间做权衡。这也是我在项目里把不确定性建模成盒式集合的根本原因：参数少、语义清楚、计算友好。

1.3 关键场景辨别在这个框架里扮演什么角色

直接求解两阶段鲁棒优化是出了名的难，因为它是一个min-max-min三层嵌套结构，理论上是一个NP-hard问题。工程上不会硬解，而是采用分解算法，比如Benders分解、列与约束生成算法。这些算法都需要在迭代过程中反复求解一个内层max-min子问题——也就是在不确定集合里寻找一个让运行成本最高的“最坏场景”。

问题在于，如果这个最坏场景是靠连续搜索在整个不确定集合里找出来的，计算量很大，而且在工程直觉上，你只需要关心那些真正可能发生的高影响场景，没必要对概率极低的极端组合做过度防御。关键场景辨别算法的思路就是：先从海量历史场景或者抽样生成的场景中，用聚类等方法筛选出一批有代表性的关键场景，用这批场景去逼近原始不确定集合中最危险的部分。这样既保留了鲁棒调度“防患于未然”的核心理念，又大幅降低了内层子问题的搜索难度，让整个求解过程在实际算例里能稳定收敛。

2. 核心原理拆解：不确定集合、关键场景辨别与鲁棒问题的求解

2.1 不确定参数的盒式集合与预算约束

在不确定性问题里，第一步是回答“不确定性长什么样”。我以光伏出力场景为例，设实际有功出力为P_pv，预测值为P_pv^0，允许的最大波动幅度为ΔP_pv^max，那么光伏出力的不确定参数可以写成：

P_pv = P_pv^0 + u_pv * ΔP_pv^max， 其中 u_pv ∈ [-1, 1]。

单独看每个参数，这是一个对称盒式区间。但这还不够——如果让所有不确定参数同时取到极端值，鲁棒优化的结果会过于保守，运行成本高到没有实际意义。真实系统里，光伏和风电的出力偏差不太可能在同一时刻都达到最小值。所以工程上引入了预算约束（Budget of Uncertainty）：

Σ_i |u_i| ≤ Γ， i ∈ Ω_uncertain。

Γ就是这个预算，它限制了不确定参数同时偏离预测值的总程度。Γ等于0时模型退化为确定性优化；Γ等于历史场景中观测到的最大总偏差时，模型退化为最保守的全盒式鲁棒。直观理解：调度员给不确定性“戴了个笼头”，不让所有坏事同时发生，这既符合物理规律，也让解决方案的保守程度可控。

2.2 两阶段鲁棒优化模型的一般形式

两阶段鲁棒模型的标准写法是：

min_x ( c^T x + max_{u∈U} min_{y∈Ω(x,u)} b^T y )。

x 是第一阶段变量，包括燃气轮机开停机状态等0-1变量；u 是不确定参数向量，落在前面定义的盒式集合U中；y 是第二阶段变量，包括各机组出力、储能充放电功率、购售电功率等连续变量；Ω(x,u) 是给定x和u后的可行域，由功率平衡、机组出力上下限、爬坡约束、储能SOC约束等构成。

这个式子的物理含义很明确：我先定第一阶段决策x，然后大自然（或者说不确定性）选择对我最不利的u，我再在已知x和u的情况下做最经济的第二阶段决策y。目标函数里的max与min是嵌套关系，反映了决策者“先决策、后观测、再响应”的顺序。

补充一点，第二阶段变量可以包含对失负荷、弃风弃光的惩罚项，这样在最坏场景下如果系统确实无法满足全部负荷，至少能在目标函数中体现代价，保证模型有解。实际代码中，我会在功率平衡约束里加入失负荷量和弃风弃光量两个松弛变量，并配一个足够大的惩罚系数。

2.3 关键场景辨别算法：先筛场景再上鲁棒

关键场景辨别的出发点很朴素的：如果我们能从海量历史场景里挑出几个“少而精”的场景，让模型在这些场景上都可行且经济，那在真实运行里大概率也不会出大问题。具体做法分三步走：

第一步，场景生成与标准化。用拉丁超立方抽样或者直接从历史数据中读取过去一整年的风光负荷出力曲线，生成大样本的场景集。注意先对数据进行标准化，消除量纲影响，否则光伏的数值范围会压过负荷的数值特征，聚类结果会失真。

第二步，聚类精简。对标准化后的场景做K-means聚类，聚类数K的选择是关键。K太小，代表性不足，最坏场景可能被“平均”掉；K太大，计算量优势又没了。我的经验是先从K=3开始，结合轮廓系数判断，一般微网的场景K取4到6比较均衡。

第三步，关键场景提取。这里有一处新手容易踩的坑：不能直接用类中心作为代表场景。因为聚类中心是类内样本的“平均画像”，它会抹平极端特征，而鲁棒优化恰恰最关心极端场景。我在项目里的做法是：先确定每个类的中心，然后在类内寻找偏离中心最大、同时总偏差水平超过某个阈值的样本，把它提取出来作为关键场景。针对每个类，至少提取三条——中心场景、偏离最大的上边界场景、偏离最大的下边界场景——共同构成关键场景集。这样既保留聚类带来的计算优势，又不丢失最坏情况的刻画力。

表1是某组测试数据下，不同关键场景数量对优化结果的影响对比，可以看出场景个数从200降到15之后，目标函数值和最终调度策略都基本稳定，但求解时间缩短了一个数量级。

仅用聚类中心的对照组目标值偏低，看起来经济性更好，但实际上丢失了边界场景，一旦真实场景出现在类边缘，方案很可能会失稳，这就是我强调要同时保留“最坏偏离样本”的原因。

2.4 求解算法：C&CG怎么把三层嵌套拆成可迭代的MILP

模型本身不能直接在Yalmip里用一条optimize命令求解，因为内层有max-min嵌套。我采用的是列与约束生成算法（C&CG）。它把原问题拆成一个主问题和一个子问题，迭代求解。

主问题（MP）是把内层max-min问题替换成一组新增约束后的MILP：

min_{x,η} c^T x + η，

约束包括第一阶段可行性约束，以及针对k个已识别场景的“可行割”：

η ≥ b^T y^k，
y^k ∈ Ω(x, u^k_hat)， k=1,2,...,K_iter。

也就是说，每找到一个新的关键场景u^k，就往主问题里加一组对应的第二阶段决策变量y^k和约束。η代表当前所有已知场景里第二阶段的最高成本下限。

子问题（SP）是给定第一阶段解x*后，求解内层max-min问题，找一个新的最坏场景u：

max_{u∈U} min_{y∈Ω(x*,u)} b^T y。

这个子问题仍然是个双层结构。如果第二阶段是线性规划，我可以用对偶理论把内层min转化为max，然后与外层max合并成一个单层max问题。合并后的SP是一个带双线性项的优化问题，不过好在这些双线性项通常是u与对偶变量的乘积，经过大M法线性化之后，可以用求解器直接求解。SP解出的最坏场景u_new_hat被加入主问题的场景集合，同时如果SP的目标函数值比当前上界还大，就更新上界，持续迭代直到上下界间隙小于设定阈值。

判断收敛条件时我用的是一次迭代里上下界的间隙，可设间隙上限为0.5%且迭代次数达到上限即停止。实测下来，C&CG相比传统Benders分解在两阶段鲁棒问题上的优势是主问题的下界收敛得更快，通常10到20轮迭代内就能稳定。

3. Matlab代码实现：整体框架与关键模块

3.1 测试系统与基础参数配置

我在项目里用的微网测试系统包含光伏、风机、蓄电池储能、微型燃气轮机和上级电网交互，共5类单元。基础参数如表2所示，这些参数会直接写入代码的参数文件里，方便批量修改。

3.2 代码目录结构与主函数流程

代码采用模块化结构，文件组织如下：

text复制two_stage_robust_microgrid/
├── main.m                    // 主脚本，控制整体流程
├── init_parameters.m         // 所有基础参数
├── data/
│   ├── pv_history.csv        // 历史光伏出力数据
│   ├── wind_history.csv      // 历史风电出力数据
│   └── load_history.csv      // 历史负荷数据
├── scenario/
│   ├── generate_scenarios.m  // 拉丁超立方抽样生成场景
│   ├── standardize.m         // 数据标准化
│   └── key_scenario_select.m // 关键场景辨别提取
├── solver/
│   ├── build_MP.m            // 构建C&CG主问题
│   ├── build_SP.m            // 构建C&CG子问题
│   └── ccg_algorithm.m       // C&CG主循环
├── utils/
│   ├── plot_results.m        // 结果绘图
│   └── save_results.m        // 结果保存
└── results/
    ├── output_x.csv          // 第一阶段决策结果
    ├── output_y.csv          // 第二阶段运行结果
    └── convergence_curve.fig // 收敛曲线

主函数流程不复杂，核心逻辑是串联数据加载、场景生成、关键场景辨别、C&CG迭代求解和结果导出。我加上时间戳计时，方便评估关键场景辨别环节相比全场景鲁棒省了多少时间。

3.3 关键场景辨别模块的代码实现

场景标准化的代码处理如下，注意一点：标准化要基于全样本的统计量，不能分批次做，否则每类的相对尺度会被破坏。

matlab复制function [scen_norm, mu, sigma] = standardize(scen_raw)
    % scen_raw: N_samples x T_scen 的原始场景矩阵
    mu = mean(scen_raw, 1);
    sigma = std(scen_raw, 0, 1);
    sigma(sigma < 1e-6) = 1; % 防止除零
    scen_norm = (scen_raw - mu) ./ sigma;
end

K-means聚类加关键场景提取的核心代码如下。聚类数我设置为4，用silhouette函数做一次质量校验，silhouette平均值低于0.5时提示用户考虑调整K值。

matlab复制function key_scen = key_scenario_select(scen_norm, K, pick_num)
    % 第1步：K-means聚类
    [idx, center] = kmeans(scen_norm, K, 'Replicates', 10, 'MaxIter', 500);

    % 第2步：在每个类内寻找偏离中心最大的样本
    key_scen = [];
    for k = 1:K
        class_data = scen_norm(idx == k, :);
        dist_vec = vecnorm(class_data - center(k, :), 2, 2);
        [~, sort_idx] = sort(dist_vec, 'descend');
        % 取偏离最大的前 pick_num 个样本，并去重
        for j = 1:pick_num
            cand = class_data(sort_idx(j), :);
            if isempty(key_scen) || min(vecnorm(key_scen - cand, 2, 2)) > 1e-3
                key_scen = [key_scen; cand];
            end
        end
    end
    % 第3步：补充类中心，保证算法的基准代表性
    key_scen = [key_scen; center];

    % 第4步：反标准化回原始物理量纲
    % 注意要用 standardize 阶段保存的 mu/sigma
end

这段代码里有一点需要特别说明：我之所以在取“偏离最大样本”之外又补充了类中心，是为了防止极端样本过于集中而失去常规运行场景的约束能力。鲁棒优化既要管住最坏情况，也不能让正常运行点完全失去经济性，两个都保留模型才会平衡。

3.4 C&CG算法主循环的代码实现

C&CG主循环里，我用Yalmip建模。第一步先构建主问题和子问题模型。

主问题模型的目标函数为：

f_mp = c_first' * x + eta

其中c_first是第一阶段成本系数，eta是标量辅助变量。每次迭代后，新识别的场景会以新约束块的形式加入主问题。先看主问题构建片段：

matlab复制function [MP, x, eta, y_cell] = build_MP(params, key_scen_set)
    yalmip('clear');
    % 第一阶段变量：燃气轮机开机状态 z_gt (24小时)
    z_gt = binvar(1, params.H, 'full');

    % 第二阶段变量池：每个已识别场景对应一组 y
    y_cell = {};   % y_cell{k} 是第 k 个场景对应的第二阶段变量结构体
    eta = sdpvar(1, 1);

    % 目标函数
    obj = sum(params.c_gt_start .* z_gt) + eta;

    constraints = [];
    % 第一阶段的启停逻辑约束 ...
    for k = 1:length(key_scen_set)
        u_now = key_scen_set{k};
        y_now = struct();
        % 第二阶段变量：燃气轮机出力、储能SOC、购售电功率、失负荷量等
        y_now.p_gt = sdpvar(1, params.H, 'full');
        y_now.p_buy = sdpvar(1, params.H, 'full');
        y_now.p_sell = sdpvar(1, params.H, 'full');
        y_now.soc = sdpvar(1, params.H + 1, 'full');
        y_now.p_loss = sdpvar(1, params.H, 'full');

        % 功率平衡约束（引入 u_now）
        constraints = [constraints, ...
            y_now.p_gt + params.p_pv_base .* (1 + u_now.pv) + ...
            params.p_wind_base .* (1 + u_now.wind) + ...
            y_now.p_buy - y_now.p_sell + y_now.p_loss == params.p_load_base + y_now.soc(2:end) - y_now.soc(1:end-1) + ...];

        % 第三阶段相关运行约束 ...
        y_cell{k} = y_now;

        % 关键：添加 η ≥ 第二段成本 的可行割
        constraints = [constraints, eta >= sum(params.c_buy .* y_now.p_buy) + ...
            sum(params.c_gt_fuel .* y_now.p_gt) + ...
            params.penalty * sum(y_now.p_loss)];
    end

    MP = optimize(constraints, obj, sdpsettings('solver', 'gurobi', 'verbose', 0));
end

子问题模型要处理的是max-min嵌套。这里我选择对第二阶段做KKT条件转换，而不是双重对偶，因为当第二阶段包含储能状态转移这类时序约束时，KKT表达在某些求解器上数值表现更稳定。代码实现如下。

matlab复制function [SP_result, u_new] = build_SP(params, x_star)
    yalmip('clear');
    % 不确定变量
    u_pv = sdpvar(1, params.H, 'full');
    u_wind = sdpvar(1, params.H, 'full');
    u_load = sdpvar(1, params.H, 'full');

    % 不确定集合约束
    constraints_U = [-1 <= u_pv <= 1, -1 <= u_wind <= 1, -1 <= u_load <= 1];
    constraints_U = [constraints_U, sum(abs(u_pv)) + sum(abs(u_wind)) + sum(abs(u_load)) <= params.Gamma];

    % 第二阶段变量
    y = struct();
    y.p_gt = sdpvar(1, params.H, 'full');
    % ... 其余变量定义省略

    % 第二阶段约束集合 Omega(x_star, u)
    constraints_Omega = [];
    % 功率平衡、储能、爬坡、购售电上限等

    % 目标函数：内层 min 部分
    obj_inner = sum(params.c_gt_fuel .* y.p_gt) + ...;

    % 通过 KKT 条件把内层 min 等价转换为互补约束
    [kkt_constraints, details] = kkt(constraints_Omega, obj_inner, y.p_gt);

    % 最终子问题：在 U 与 KKT 条件下最大化内层目标
    obj_sp = -obj_inner;  % 因为kkt求解的是min问题的稳定点，再对外取max
    optimize([constraints_U, kkt_constraints], -obj_sp, sdpsettings('solver', 'gurobi', 'verbose', 0));
end

严格来说，KKT条件转换后的互补约束是非线性非凸的，求解器不一定直接处理。工程上更稳的做法是使用对偶重写，将内层min问题转化为对偶变量下的max问题，得到的SP是一个带双线性项的单层优化问题。双线性项形如u * λ，由于u和λ都是有界变量，可以引入辅助变量并用大M法精确线性化。这块代码比较长，但逻辑是固定的：写出对偶约束、对偶目标，然后把双线性项替换为辅助变量与分段线性约束。

3.5 C&CG迭代主循环与收敛判据

主循环部分逻辑如下：

matlab复制function [result, history] = ccg_algorithm(params, scen_init)
    LB = -Inf; UB = Inf;
    history.ub = []; history.lb = [];
    key_scen_set = scen_init; % 初始关键场景

    gap_tol = params.gap_tol;   % 比如 0.005
    max_iter = params.max_iter; % 比如 30

    for iter = 1:max_iter
        % 1. 求解主问题
        [x_opt, obj_mp, y_set] = solve_MP(params, key_scen_set);
        LB = max(LB, obj_mp);

        % 2. 求解子问题，得到最坏场景
        [u_new, obj_sp] = solve_SP(params, x_opt);
        UB = min(UB, obj_sp);

        % 3. 判定收敛
        gap = (UB - LB) / max(1e-6, abs(UB));
        history.ub(end+1) = UB;
        history.lb(end+1) = LB;

        if gap <= gap_tol
            break;
        end

        % 4. 把新场景加入主问题集合
        key_scen_set{end+1} = u_new;
    end

    result.x = x_opt;
    result.UB = UB;
    result.LB = LB;
    result.iter = iter;
    result.key_scen_num = length(key_scen_set);
end

这里我遇到过一个很典型的现象：UB和LB初始值设置不当会导致首轮gap异常。初始UB应该设为Inf，初始LB设为-Inf，第一轮主问题解出来后LB更新为有限值，第一轮子问题解出来后UB才更新。如果你的代码首轮gap输出NaN，先检查这两个初始值。

4. 运行结果与调度行为分析

4.1 三种典型天气下的调度策略

我用三条典型日曲线做了测试：晴朗日（光伏高出力）、阴天（光伏低出力、风电中高出力）、极端晚峰（负荷晚高峰叠加风电低出力）。

晴朗日的调度策略比较有代表性。第一阶段燃气轮机开启时段明显减少，储能SOC保持在中低水平，为午间光伏大发预留充电空间。到了午后，光伏出力进入高位，储能开始充电，燃气轮机停机。晚间负荷上升，光伏出力归零，储能放电、燃气轮机启动、适时购电。鲁棒解和确定性解在这个场景下差距不大，因为预测误差相对小，最坏场景带来的额外成本有限。

阴天场景是鲁棒优化的价值体现。模型预测光伏出力只有晴天的四成，同时不确定性预算发挥作用——鲁棒解中燃气轮机的启停次数比确定性解多一次，储能从早上开始就保持较高SOC，为可能的光伏出力不及预期预留充足调节空间。购电功率上限在晚峰时段也更早触发，说明模型在最坏情况下优先保证负荷供应。

极端晚峰场景最考验方案质量。在这个场景下，确定性方案晚峰时段会出现约45kW的失负荷量，而鲁棒方案通过提前一小时启动燃气轮机、储能放电功率上限预留10%裕量，把失负荷降为0。代价是总运行成本比确定性方案提高了约11%。这个成本上升可以理解为“买保险”的保费。

4.2 迭代收敛行为与关键场景自更新过程

图1展示了某次测试的C&CG收敛曲线，UB与LB在8轮迭代内从初值逼近到间隙小于0.5%。值得关注的是，第2轮迭代识别出的最坏场景不是光伏最低的一小时，而是“光伏连续三小时低出力+晚负荷抬升”的组合场景，这说明关键场景辨别在算法迭代中会自动向系统最脆弱的时段偏移。单独看任一小时的偏差，这个场景并不算极端，但它触发了储能SOC的全局约束，对调度方案的影响远大于单小时极端场景。

这个现象引出一个实操总结：两阶段鲁棒问题里真正危险的不确定场景往往是“时间组合型”的，例如连续多小时的出力偏低叠加晚高峰负荷上升。关键场景辨别算法通过聚类和偏离度筛选，天然保留了这种组合特征，而全场景抽样方法反而容易在高维抽样中错过这类组合场景的小概率区间。

4.3 关键场景数量与求解质量的权衡

我把关键场景数量从3逐步增加到24，观察目标函数值和求解时长的变化。结果显示，当关键场景数从3增加到6时，目标函数值有明显跳升（因为最坏场景被纳入了考虑），从6增加到18时目标值缓慢上升，超过18后基本趋于平稳。求解时长则几乎随场景数线性增长。

从工程角度，我建议关键场景数量控制在总场景数的5%到10%之间，并通过灵敏度测试判断是否已经覆盖最坏场景组合。当关键场景数量超过15后，目标值变化小于0.3%时，就没必要再继续增加了。

5. 常见问题与调试实录

5.1 内层子问题对偶转换失败或不可行

这是我自己踩过最多次的坑。用对偶方法重写第二阶段问题时，如果第二阶段原始问题不可行，对偶问题也无解，整个SP求解会报 infeasible。排查思路很明确：先单独固定一组u值，解一个纯第二阶段问题，确认原始问题是否可行；再把对偶约束逐条检查，确认是否漏写了某一组对偶约束。另一个高频原因是第二阶段里存在等式约束，对偶变量没有符号限制，但代码里写成了非负约束，导致对偶可行域被错误缩小。

处理办法：给第二阶段功率平衡约束加非负松弛变量，即失负荷量和弃风弃光量，并设置适当大的惩罚系数。这样即使u取到某类特殊组合，原始问题也一定有可行解，对偶问题随之可解。同时，第二阶段约束里避免使用严格不等式，所有上下限都用≤、≥，否则对偶推导会出问题。

5.2 Yalmip报错“Infeasible problem”但检查不出原因

如果MP或SP单独求解时报不可行，而约束看起来都很合理，我建议按以下顺序排查：

1. 检查第二阶段变量维度是否与调度周期一致。储能SOC变量维度经常出错，SOC_0的初始化位置和SOC_H+1的终值约束要仔细核对。
2. 检查功率平衡约束中换流器损耗的处理。交流微网里DC/AC换流器有损耗，如果忽略会导致功率无法平衡。
3. 检查极端场景下所有可调资源的总上限是否小于负荷峰值。比如光伏出力为0、风电出力为0、燃气轮机已达到上限、购电也达到上限情况下，系统总出力是否还能满足负荷。如果不满足，模型必然不可行，这不是代码问题，而是系统容量本身就不足，需要增加松弛变量或调整设备配置。

5.3 C&CG迭代次数过多，收敛过慢

场景预算Γ设置过小时，内层搜索空间变小，可能很快收敛；但当Γ较大且第二阶段包含储能时序约束时，迭代次数会显著增加。我的经验性调整参数如下：

另外，Gurobi和Cplex对MILP的求解效率有明显差异，同一模型在Gurobi上比Cplex快约30%。如果碰到求解时间离谱，先换个求解器试试。

5.4 求解器许可证与安装问题

这套代码里求解MILP和MIQP的求解器不是Matlab自带的。我主要用Gurobi和Yalmip搭配。如果只装了Cplex，会遇到许可证激活失败、找不到求解器的问题。建议直接安装Gurobi并申请学术许可证，注册edu邮箱即可，配置时执行gurobi_setup。另外，求解器安装完必须在Matlab里运行yalmiptest，等输出列表中看到Gurobi对应的MPC、MILP等条目显示可用，再跑主程序。

5.5 结果画图与调度曲线可视化

结果展示建议至少画三张图：各机组24小时出力堆叠图、储能SOC曲线、C&CG收敛曲线。堆叠图用area函数画，能直观看到每个时段各类电源的贡献比例。SOC曲线注意要单独画在右侧纵轴，否则量纲差异会让曲线看不清。收敛曲线用semilogy画上下界间隙的对数曲线，比线性坐标更直观地展示收敛速度。

6. 从复现到扩展：这个方案还能怎么改

两阶段鲁棒微网调度这个课题做好了，扩展空间非常大。最简单的一块是给储能加上容量衰减成本，把电池的循环老化折算成单位充放电成本，这样调度策略就不会把储能当成免费资源随意使用。另一个方向是引入需求响应资源，把可中断负荷作为第三道防线，参与最坏场景下的功率平衡。再往深走，可以在第二阶段加入电动汽车集群的充放电调度，考察交通行为与电力系统耦合下的鲁棒调度问题。

如果是从零开始复现这套代码，我建议严格按三个阶段推进。第一阶段先做确定性版本，确保微网模型本身没有错误。第二阶段加入关键场景辨别，用15个关键场景替代全场景集合，观察目标值变化。第三阶段再引入C&CG迭代，实现完整的鲁棒求解。跨过第一阶段的代码问题直接写鲁棒版，出了问题根本分辨不清是物理模型的问题还是算法逻辑的问题。我自己第一次复现时就在功率平衡约束里漏了储能SOC的时序关系，结果SP里的最坏场景始终集中在同一时段，收敛曲线一直震荡，排查了一个下午才定位到根因。

在大大小小的微网调度项目里转过一圈之后，我最大的体会是：鲁棒优化解决的不是“预测得准不准”的问题，而是“预测不准时系统会不会崩”的问题。关键场景辨别算法在这个框架里就像一道筛子，把海量不确定性的细沙滤掉，只留下真正决定调度方案成败的那几块石头。把这套主逻辑打通之后，换成其他不确定性来源、其他微网拓扑、甚至其他两阶段决策问题，剩下的都只是改参数换约束的体力活。